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ABSTRACT 


Six  vertical  discretization  schemes  are  compared  for  a 
linear,  baroclinic,  vorticity-divergence  equation  model. 
Variables  are  defined  for  two  staggered  and  one  unstaggered 
grids.  -A  finite  difference  and  a  Galerkin  finite  element 
approximation  are  formulated  for  each  of  the  three  grids. 
The  models  are  tested  in  three  experiments.   The  largest 
difference  between  the  grids  occurs  in  the  mid-atmosphere 
diabatic  heat  source  experiment  and  the  unstaggered  grid 
provides  the  best  results.   For  the  staggered  grids,  the 
results  of  the  finite  element  models  are  not  more  accurate 
than  the  corresponding  finite  difference  results.   Oscilla- 
tions occur  in  the  temperature  profiles  near  the  surface  for 
both  staggered  finite  element  models. 
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I.   INTRODUCTION 

The  finite  elements  method  (FEM)  has  been  used  in  several 
atmospheric  prediction  models  (Staniforth  and  Mitchell,  1977; 
Williams  and  Schoenstadt,  1980,  Beland,  et  al.,  1983).   The 
FEM  is  a  special  case  of  the  Galerkin  procedure  in  which  the 
dependent  variables  are  approximated  by  a  finite  sum  of 
spatially  varying  basis  functions  with  time  dependent  coeffi- 
cients.  The  FEM  basis  functions  are  low  order  polynomials 
which  are  zero  except  in  a  localized  region.   The  Galerkin 
procedure  produces  a  set  of  coupled  ordinary  differential 
equations  which  are  solved  by  using  finite  differences  in 
time  (Haltiner  and  Williams,  1980) . 

FEM  models  are  potentially  more  accurate  than  finite 
difference  method  (FDM)  models.   Winninghoff  (1968),  Arakawa 
and  Lamb  (1977),  and  Schoenstadt  (1980)  have  demonstrated 
the  advantages  of  spatial  staggering  of  predictive  variables 
in  finite  difference  approximations  of  shallow  water  equa- 
tions.  The  results  of  Williams  and  Schoenstadt  (1980)  indi- 
cated that  FEM  models  should  either  use  staggered  nodal  points 
in  the  momentum  equations  or  unstaggered  nodal  points  in  the 
vorticity-divergence  equations. 

The  purpose  of  this  research  is  to  compare  six  linear, 
baroclinic,  vorticity-divergence  equation  models.   A  finite 
difference  and  a  finite  element  model  are  written  for  each  of 
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the  three  vertical  grids  depicted  in  Figure  1.   Grid  A,  which 
was  originally  developed  by  Lorenz  (1960)  for  energy  conser- 
vation, is  a  widely  used  grid  for  finite  difference  models. 
However ,~ Tokioka  (1978)  has  shown  that  this  grid  has  a 
computational  mode  in  the  temperature  field,  and  Arakawa  (19  84) 
has  found  a  false  small  scale  baroclinic  instability  with  this 
grid.   Tokioka  (1978)  also  found  that  grid  C  has  computational 
modes  in  all  fields.   Beland,  et  al.,  (1983),  analyzed  the 
finite  element  formulation  for  grid  C  and  found  noise  gener- 
ated by  certain  forms  of  friction.   Grid  B  was  introduced  by 
Charney  and  Phillips  (1953)  for  a  finite  difference  quasi- 
geostrophic  model.   Tokioka's  (1978)  analysis  has  shown  that 
this  grid  has  no  computational  modes.   It  is  hoped  that  the 
finite  element  formulation  for  this  grid  will  have  high 
accuracy  without  the  problems  the  other  grids  have  had  with 
computational  modes. 

Each  of  the  models  used  in  this  thesis  are  derived  from 
the  governing  equations  described  in  Chapter  II.   The  models 
are  one-wave  spectral  in  the  horizontal  and  have  a  fixed  lid 
at  the  top  of  the  atmosphere.   Diabatic  heating  and  mountain 
topography  can  be  included  in  each  model  experiment. 

The  results  of  three  experiments  with  analytic  initial 
conditions  will  be  examined.   As  described  in  Chapter  III,  the 
experiments  are  1)  an  initial  perturbation  in  the  meridional 
flow,  2)  flow  over  mountains,  and  3)  flow  with  a  mid-atmosphere 
diabatic  heat  source.   Each  experiment  will  be  presented  with 
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Figure  1.   Three  vertical  grids  to  be  compared  for  N- 
layer  models.   Perturbation  variables  are 
defined  at  either  unstaggered  levels  (solid 
lines  at  height  Z')  or  staggered  levels 
(dashed  lines  at  height  Z)  .  z,    is  vorticity, 
D  is  divergence,  T  is  potential  temperature, 
T  is  basic  state  potential  temperature, 
<j)  is  geopotential ,  <j>s  is  surface  geopotential, 
u  is  east-west  velocity,  u  is  basic  state 
east-west  velocity,  v  is  north-south  velocity, 
Z  is  vertical  velocity,  Q  is  dibatic  heating, 
and  MTS  is  forced  vertical  velocity. 


14 


six  and  sixty  vertical  layers.   The  results  of  the  60-layer 
models  will  show  if  the  models  are  converging  to  the  same 
solution.   For  each  model,  comparison  of  the  six  and  60- 
layer  results  indicate  how  well  the  lower  resolution  model 
approximates  the  convergent  solution.   The  purpose  of  these 
experiments  is  to  identify  any  problem  areas  in  the  models 
and  to  demonstrate  the  general  characteristics  of  each  model 
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II.   MODEL  DESCRIPTIONS 

A.  MODEL  FEATURES 

The  six  numerical  models  are  developed  with  several  fea- 
tures which  make  the  models  easy  to  modify  for  a  wide  range 
of  experiments.   The  depth  of  the  atmosphere,  number  of 
levels  and  the  depth  of  each  layer  are  variable  in  each  model. 
Diabatic  heating  and  forced  vertical  velocity  due  to  mountain 
topography  are  included  in  the  governing  equations.   The  user 
can  prescribe  heating  or  forced  vertical  velocity  functions, 
or  run  the  model  with  these  terms  defined  as  zero.   The  models 
are  written  in  modular  structure  using  FORTRAN  '77.   There  is 
parallel  construction  between  models.   The  subroutines  used 
in  one  model  are  very  similar  to  those  used  in  the  other  five 
models.   The  models  run  very  quickly  on  an  IBM-3033  mainframe. 
A  96-hour  forecast  for  a  12-layer  FEM  uses  less  than  five 
seconds  of  computer  processing  time. 

B.  govf:rning  EQUATIONS 

Each  model  approximates  the  same  set  of  governing  equa- 
tions.  The  vorticity  (2-1),  divergence  (2-2),  surface  geo- 
potential  (2-3)  equations  and  the  first  law  of  thermodynamics 
(2-4)  are  the  prognostic  equations  for  the  forecast  variables 
vorticity,  divergence,  surface  geopotential  and  potential 
temperature.   The  surface  geopotential  equation  is  the  lower 
boundary  condition  on  the  vertical  velocity.   The  vertical 
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coordinate  Z  =  -ln(p/pQ)  is  used,  but  the  non-Boussinesq 
terms  involving  e  are  replaced  by  one.  The  prognostic 
equations  in  the  coordinates  x,y,Z,t  are 


d£    /   .exr-v    „     3Z  3v    3Z  3u      n  ,„  ,, 

=±    +   (£+f)D  +  $V  +  ~—  -T ^—  -r—   =    0,  2-1 

dt  3x  3x    3y  3y 


dD    ,3u. 2    ,9v,2    .  3u3v    3Z  3u 
dt    l3x;     l3y;     z  3y  3x   9x  3Z 

+  ||  ||  +  3u  -  fc  +  V2(J)   =   0  f  (2-2) 

d$_ 

=   MTS  ,  (2-3) 


dt 
dT 


dt   =   Q  .  (2-4) 


Here, 


€,      is  the  vertical  component  of  vorticity,  £  =  t— ■  -  -r— , 

D   is  the  horizontal  divergence,  D  =  tt—  +  t—, 

3  3x    3y' 

$  is  the  geopotential ,  $  =  gZ, 

(J)  is  the  surface  geopotential, 

T  is  the  potential  temperature, 

u  is  the  x-component  of  velocity, 

v  is  the  y-component  of  velocity, 

Z  is  the  vertical  velocity, 

Q  is  the  diabatic  heating  per  unit  time  per  unit  mass, 

MTS  is  the  forced  vertical  velocity  due  to  flow  over 
mountain  topography,  discussed  in  Chapter  III, 
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f   is  the  Coriolis  parameter, 

6  1S  dy' 

d     9       9       3     19 
dt    9t      9x      9y      9Z 

5  "  ■ 
V   is  the  horizontal  Laplacian  operator, 


The  prognostic  equations  are  linearized  by  expanding  the 
variables  into  the  following  mean  state  and  perturbations: 

u(x,Z,t)   =   u(Z)  +  u'(x,Z,t)  ,  (2-5) 

v(x,Z,t)   =   v' (x,Z,t)  ,  (2-6) 


Z(x,Z,t)   =   Z(x,Z,t)  ,  (2-7) 


T(x,y,Z,t)   =   T(y,Z)  +T'(x,Z,t)  , 


<j>(x,y,Z,t)   =   4>(y,Z)  +  $l(x,Z,t)     ,  (2-9) 


<j>s(x,y,t)   =   «J>s(y)  +  <J>^(x,t: 


C(x,Z,t)   =   £'(x,Z,t)  , 


D(x,Z,t)   =   D' (xf Z,t)  , 


Q(x,Z,t)   =   Q'(x,Z,t)  , 


MTS (x,t)   =   MTS' (x,t) 
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3v' 
The  vorticity  and  divergence  in  this  system  are  c,   =    ~ — , 

3u' 
D  =  x •   The  linearized  forecast  equations  are 

a   X 


If    ="  -fD'   "  u  If  '   Bv'  '  (2"15> 

3x 


3T 


-  3T'     ,  3T    I  3T 


3t         "  3x        3y      3Z 


Q'  ,  (2-17) 


3  4)'  34)'       3$      _. 

=   -  u  ,r-2.  -  v'  ■*— -  -  RTZ  +  MTS'  ,  (2-18) 


3t  3x        3y 

where  R  is  the  gas  constant  for  air. 

The  diagnostic  variables,  u'  ,  v'  ,  Z',  <p '  ,  are  calculated 
from  the  forecast  variables  using  the  definitions  of  divergence, 
vorticity,  the  hydrostatic  equation  and  the  continuity  equa- 
tion.  The  relations  are 


3u' 

3x 


3v' 
3x 


=   D'  ,  (2-19) 


=  V     ,  (2-20) 


34>' 
3x 


RT'  ,  (2-21 


D'  +  ||1   =   0  .  (2-22) 


The  use  of  primes  to  denote  perturbation  quantities  will  be 
discontinued.   All  quantities  used  in  the  remainder  of 
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the  paper  will  be  perturbation  quantities  unless  otherwise 
noted. 

The  mean  state  is  assumed  to  be  in  hydrostatic  and  geo- 
strophic  balance.   The  term  3T/3y  in  the  first  law  of  thermo- 
dynamics can  be  evaluated  by  taking  3/9y  of  the  hydrostatic 
equation  and  substituting  for  d<$>/dy   from  the  geostrophic 
relation,  3<J>/3y  =  -fu.   Thus, 


21   =   _  f  9u 

3y        R  9Z  *  (Z    ZJ; 


Geostrophic  balance  of  the  mean  state  at  the  surface  implies 

=   -  fu  .   .  (2-24) 


9y  sfc 

The  expressions  (2-23)  and  (2-24)  are  substituted  into  equa- 
tions (2-17)  and  (2-18),  respectively. 

A  singlewave  spectral  representation  is  used  in  the 
x-direction,  with  wave  number  y  =  2tt/L,  where  L  is  the  wave- 
length in  the  x-direction.   The  perturbation  quantities  have 
the  form 

C'(x,Z,t)   =   A,  (Z,t)cos  yx  +  A2(Z,t)sin  yx  ,  (2-25) 

D'(x,Z,t)   =   D1(Z,t)cos  yx  +  D2(Z,t)sin  yx  ,  (2-26) 


T'(x,Z,t)   =   T-,(Z,t)cos  yx  +  T2(Z,t)sin  yx  ,  (2-27) 
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<J>'(x,t)   =   S1(t)cos  yx  +  S2(t)sin  yx  ,  (2-28) 


u' (x,Z,t)   =   U1(Z,t)cos  yx  +  U2(Z,t)sin  yx  ,  (2-29) 


v'(x;Z,t)   =   V1(Z,t)cos  yx  +  V2(Z#t)sin  yx  ,  (2-30) 


Z'(x,-2,t)   =  W1(Z,t)cos  yx  +  W2(Z,t)sin  yx  ,  (2-31) 


<J>'  (x,Z,t)   =   H,(Z,t)cos  yx  +  H2(Z,t)sin  yx  ,  (2-32) 


Q'  (x,Z,t)   =   (^(Z/tJcos  yx  +  Q2(Z,t)sin  yx  ,  (2-33) 


MTS 


1 (x,t)   =   MTS1(t)cos  yx  +  MTS2(t)sin  yx  .  (2-34) 


The  relations  (2-25) -  (2-34)  are  substituted  into  equations 
(2-15) - (2-22) .   The  prognostic  and  diagnostic  equations  are 
separated  into  equations  for  the  cosine  and  sine  terms.   The 
resultant  prognostic  equations  are 


dA 

j^-     =      -fD1  -  uyA2  -  3V1  ,  (2-35) 


dA 

j^l      =   -fo2  +  uyAx  -  6V2  ,  (2-36) 


at1    =    fAi  "  U^D2  "  6Ui  "  »  a£  W2  +  ^Hl  '  (2"37) 
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8D2  _  du       2 

y^-   =   fA2  +  uuD1  -  6U2  +  y  ^  Wx  +  U  H2  ,  (2-38) 


1       —       f  du        9T 
1t~     ="  _UUT2  +  R  dZ  Vl  "  TJz  Wl  +  °1  '  <2"39> 


gm  

TtT  -   ^Tl  +  I  aT  V2  -  H  W2  +  Q2  '  -         (2"40> 


3S1       -       - 

g^   =   -uyS2  +  fuV,  -  RTW,  +  MTSX  ,  (2-41) 

=   uyS,  +  fuV0  -  RTWn  +  MTS~  .  (2-42) 


3t        "1       2       2       2 


The  resultant  diagnostic  equations  for  u  and  v  are 


D2 
U,   =   "  —  ,  (2-43) 

1        u 


(2-44) 


(2-45) 


2-46 


Geopotential  values  above  the  surface  are  obtained  by  inte- 
grating the  hydrostatic  equation  from  the  surface  (Z  =  Z~)  to 
height  Z: 


Z 
H    =   R    /   T, (Z,t)dZ  +  S,  (2-47) 

Z0 
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U2 

= 

Dl 

Vl 

= 

A2 

V2 

= 

z 

H2   =   R    /   T2(Z,t)dZ  +  S2  (2-48) 

Z0 


The  vertical  velocity  is  calculated  by  integrating  the  con- 
tinuity equation  from  the  top  of  the  atmosphere  (Z  =  Z  )  down 
to  height  Z.   The  upper  boundary  condition,  Z  =  0  at  Z  =  Z  , 
is  used.   This  boundary  condition  is  not  exact,  but  some  form 
of  it  is  used  in  most  numerical  models.   The  diagnostic 
equations  for  the  vertical  velocity  are 


ZT 
W,   =    /    D, (Z,t)dZ  ,  (2-49) 

Z 


ZT 

W9   =    /    D9(Z,t)dZ  .  (2-50) 

-    Z      A 


Equations  (2-35) - (2-50)  are  the  prognostic  and  diagnostic 
equations  that  govern  all  six  numerical  models.   Using  the  given 
basic  state  and  the  one-wave  spectral  perturbation  quantities, 
the  governing  equations  reduce  to  functions  of  Z  and  t.   The 
models  are  effectively  one-dimensional. 

To  display  the  results  of  each  model,  the  sine  and  cosine 
amplitudes  of  each  variable  are  combined  to  determine  the 
amplitude  and  phase  of  a  single  cosine  wave  in  the  x-direction. 
A  typical  variable  has  the  form 

Y(x,Z,t)   =   A(Z,t)cos(ux-6)  ,  (2-51) 
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where  the  amplitude  is  A(Z,t)  and  the  phase  is  6(Z,t).   The 
amplitude  and  phase  are  calculated  at  each  level  for  all 
variables. 

C.   TIME  DIFFERENCING 

Two  forward  time  steps  are  taken  to  start  each  model  and 
then  leapfrog  time  differencing  is  used.   The  leapfrog  scheme 
is  employed  because  of  its  ease  to  code.   A  Robert  filter  is 
used  to  reduce  the  amplitude  of  the  computational  mode  gener- 
ated by  the  leapfrog  time  differencing.   The  filter  is  dis- 
cussed by  Haltiner  and  Williams,  1980.   For  a  prognostic 

variable  F,  calculate  F   , ,  the  average  value  of  F  at  time 

n-1  3 

step  (n-1) At,  using  equation  (2-52).   Using  the  unaveraged 

3F 
values  at  time  step  nAt,  compute  the  tendency  (^-r-)   from  its 

dt  n 
predictive  equation.   The  predicted  value  at  time  step 

(n+l)At  is  calculated  using  equation  (2-53). 

F„  i   =   F   ,  +  Y(F„  -2F^  ,  +F^   )  ,  (2-52) 

n-1      n-l      n    n-1    n-2 

where  y    is  a  weighting  function. 


F  _,,   =   F   .  +  2At(.|&  (2-53 

n+1      n-1        9t  n 


In  all  thesis  experiments,  y   =    0.05  is  used.   The  time  step 
for  each  experiment  is  calculated  in  the  model  by  requiring, 
for  computational  stability, 
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vAt   -  j   ,  (2-54) 


where  v  =  yc,  and  c  is  the  typical  phase  speed  of  an  external 
gravity  wave. 

D.   VERTICAL  GRIDS 

Each  of  the  models  uses  one  of  three  vertical  grids.   The 
three  ways  of  distributing  the  variables  over  discrete  levels 
are  depicted  in  Figure  1.   The  staggered  levels  are  represented 
by  the  dashed  lines  in  Figure  1.   Notice  that  the  heights  at 
which  the  variables  are  defined  change  between  the  three 
grids.   The  notation  used  in  this  paper  to  denote  the  staggered 
and  unstaggered  levels  is  consistent  with  the  conventions  used 
in  the  coded  models.   The  height  of  the  unstaggered  levels 
is  denoted  as  Z'.   The  height  of  the  staggered  levels  is 
denoted  as  Z.   In  the  models,  both  Z'  and  Z,  are  defined  to 
be  the  surface  of  the  earth.   It  is  assumed  that  the  staggered 

level  Z-  is  exactly  in  the  middle  of  the  layer  between  Z!  , 

i  J  l-l 

and  Z!.   This  distinction  is  important  because  the  models  can 
have  layers  with  unequal  depth.   Thus,  the  height  of  the 
staggered  levels  is  defined  relative  to  the  height  of  the 
unstaggered  levels. 

A  finite  difference  model  is  written  for  each  of  the  grid 
structures.   The  models  are  denoted  as  DFM-A,  FDM-B  and 
FDM-C.   Similarly,  finite  element  models  using  the  three  grids 
are  indicated  by  FEM-A,  FEM-B  and  FEM-C. 
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E.  FINITE  DIFFERENCE  MODELS 

The  only  differences  in  the  equations  between  the  three 
FDM  models  are  the  approximations  of  terms  involving  du/dZ 
and  8T/9Z  in  the  prognostic  equations  and  the  approximations 
of  the  integral  in  the  diagnostic  geopotential  equation. 
Centered  difference  approximations  are  used,  except  at  the 
boundaries  where  one-sided  differences  are  employed.   The 
finite  difference  approximations  used  in  the  prognostic  equa- 
tions are  listed  in  Appendix  A. 

F.  FINITE  ELEMENT  MODELS 
1.   FEM-C 

The  unstaggered  FEM  model  is  the  simplest  of  the  three 
FEM  models.   Each  of  the  dependent  variables  is  expanded  into 
a  finite  series  in  terms  of  the  eigenfunctions  cJ).(Z).   The 
eigenf unctions  for  this  model  are  depicted  in  Figure  2.   The 
eigenfunction  expansion  for  a  typical  term  is 


N+l 
A1(Z,t)   =    I      Aj(t)<t>.  (Z)  .  (2-55) 

j  =  l 


The  finite  element  approximations  for  the  vorticity, 
divergence  and  thermodynamic  equations  are  derived  by  substi- 
tuting the  eigenfunction  expansion  for  each  dependent  variable 
into  equations  (2-35) -  (2-40)  .   Each  equation  is  multiplied 
by  <J>-  (Z)  and  integrated  with  respect  to  Z  from  the  bottom  to 
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Z' ■  t  i  i  i  i  i  i  i 

N  +  1 


ZN 


Z2 


l\         I  I   I   I   I   II   I   I   I   I   I   I   I 


Figure  2.   Eigenf unctions  defined  at  unstaggered  levels 
(solid  lines)  for  grid  C.   Dashed  lines 
represent  staggered  levels. 


the  top  of  the  atmosphere.   Each  term  in  the  equations  is  the 
finite  sum  of  separate  integrals.   Only  the  integrals  of 
overlapping  eigenfunctions  are  nonzero.   The  resultant  equa- 
tions, listed  in  Appendix  B,  are  matrix  equations.   For  an 
N-layer  model,  the  variables  are  calculated  at  N+1  discrete 
levels  and  the  variable  can  be  represented  as  a  column  vector 
of  length  N+1.   Thus,  the  forecast  equations  for  vorticity, 
divergence  and  potential  temperature  become  vector  equations 
which  contain  four  (N+1)  *  (N+1)  matrices  multiplying  appro- 
priate vectors. 


27 


The  mass  matrix  for  this  model  is  defined  by 


i+1      ZT 

M  =    I  j        <f>.  (Z)<(>.  (Z)dZ,   for  i  =  1,...,N+1  .     (2-56) 

J-i-i  z0    3 


The  matrix  N  is  defined  for  terms  multiplied  by  u, 


2 

i+1    i+1  T 


N   = 

j 


I  I         U.    /    <}>.  (Z)<f>.(Z)<j>.  (Z)dZ, 

=i-l  k=i-l   3  ZQ     D    k    i 

for  i  =  1, . . . ,N+1  . 
The  matrix  P  is  defined  for  terms  multiplied  by  du/dZ, 


i+1   . i+1        ZT  d0 . 
1=1  I        u.    /   -^  4)k(Z)<D.(Z)dZf 

j=i-l  k=i-l   D  ZQ     dZ    K     1 

for  i  =  1 , . . . ,N+1  . 


The  matrix  R  is  defined  for  terms  multiplied  by  9T/3Z, 


2 

i+1    i+1         T  d<|> . 

R  -    I  I       t.   /   g=l  <(>k(z)(Di(z)dz, 

j=i-l  k=i-i  :  zQ    dZ       x 

for  i  =  1 , . . . ,N+1  . 


The  method  to  evaluate  the  terms  of  matrices  (2-56) -  (2-59) 
and  the  resultant  formulas  for  the  matrix  elements  are  listed 
in  Appendix  C. 
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The  vorticity,  divergence  and  thermodynamic  equations, 
written  in  matrix  (indicated  by  a  double  underline)  and 
vector  (indicated  by  a  single  underline)  form  are 


M~(A1)   =   M-  (-f  -Dl  -  8-V1)  -  y-N-A2  ,  (2-60) 


=  dt 


M-§r-(.A2)   =   M-  (-f -D2  -  3'V2)  +  yN-Al  ,  (2-61) 


J3 

M-|^-(D1)   =   M-  (f -Al  -  3UD  -  U-N-D2  -  y  P  •  W2 

+  y2-M-Hl  ,  (2-62) 


M~(D2)   =   M-  (f -A2  -  6U2 )  +  y-N-Dl  +  yP-W2 

+  y2-M-H2  ,  (2-63) 


M~(T1)   =   -yN-T2  +  (|)  -P-Vl  -  R-Wl  +  M-Ql  ,  (2-64) 


M-~(T2)   =   yN-Tl  +  (|)  -P-V2  -  R-W2  +  M-Q2  .  (2-65) 


Equations  (2-60) -  (2-65)  are  simplified  by  multiplying 
each  equation  by  M   and  applying  the  Robert  filter.   The 
matrices  M   -N,  M   -P,  M   *R  are  constants.   They  are  con- 
structed in  the  initialization  subroutine  and  stored  for  use 
in  the  forecast  subroutine.   The  matrices  are  multiplied  by 
the  appropriate  vectors  with  values  for  time  level  nAt.   The 
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resultant  forecast  equations  are  vector  equations  and  the  fore- 
cast value  for  the  i-th  vertical  level  is  the  sum  of  values 
in  the  i-th  location  of  each  vector  in  the  equation.   The 
prognostic  equations  for  the  vorticity,  divergence  and  poten- 
tial temperature  vectors  are 


Al  L,  *  ■   A1BAR   ,  +  2At(-f *D1  -  3-V1  -  u-M~  -N'A2) 
— n+1     n-1  —   M  —   H  =====  —  n 


A2  .n   =   A2BAR   ,  +  2 At (-f • D2  -  6 • V2  + y -M*1 -N- Al )   ,      (2-67) 
— n+J.      n— x  —    —        —  —  j} 


Dl   -   =   DlBARn    +  2At(f -Al  -  g'Ul  -  y'M  1 • N - D2 


-  yM~1-P-W2  +  y2Hl)   , 

=====  n 


D2   ,   =   D2BAR   ,  +  2 At ( f • A2  -  6  * U2  + y • M~  -N-Dl 
— n+1      n-1  —     —   K  =====  — 

+  y'M   *P'W1  +  p2H2_) 


Tl  ,,   =   T1BAR   .  +  2At(-yM~  -N-T2 
— n+1      n-1         M  =====  — 

+  (& *M"  -P-Vl  -  M~  -R-Wl  +  Ql) 
R  ======  —   ======  —   —  n 


T2       ,       =      T2BAR      ,     +    2At(yM       •  N  •  Tl 
— n+1  n-1  M    =====  — 

+     (£)-M~1-P-V2    -    M~1-R-W2    +    Q2) 
R      =zzr=zrr   ====  n 


where  the  subscripts  n+1,  n,  n-1,  refer  to  the  values  of  the 
vectors  at  time  step  (n+1) At,  nAt,  and  (n-1) At,  respectively 
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The  surface  geopotential  and  the  diagnostic  variables 
are  calculated  using  the  same  equations  that  are  used  in 
the  unstaggered  model  FDM-C. 
2.   FEM-A 

The  FEM-A  model  defines  vertical  velocity  at  the  un- 
staggered levels  in  terms  of  the  eigenf unctions  \p  .  (Z)  .   The 
other  variables  are  defined  at  the  staggered  levels  in  terms 
of  the  eigenf  unctions  <j> .  (Z)  .   The  eigenf  unctions  for  this 
model  are  depicted  in  Figure  3. 


z   ''''»'' 

N+1 


ZN*T 


2N 


ZN 


Z3 


Z2 


Z2 


l\  I   I   I   I   I   I   I   I   I   I   I   I   I   I 


Figure  3.   Eigenf unctions  for  grids  A  and  B.  Egienf unc- 
tions i|>(Z)  are  defined  for  the  unstaggered 
levels  (solid  lines  at  height  Z')  and  tj>(Z) 
are  defined  for  the  staggered  levels 
(dashed  lines  at  height  Z) . 
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The  finite  element  approximations  for  the  vorticity, 
divergence  and  thermodynamic  equations  are  derived  using  the 
method  described  for  the  FEM-C  model.   The  resultant  Galerkin 
formulation  of  the  equations,  listed  in  Appendix  D,  are  matrix 
equations.   The  mass  matrix,  M,  and  the  matrix,  N,  for  terms 
of  the  form  u  multiplied  by  variables  defined  at  the  nodal 
points  Of  <j>  .  (Z)  have  the  same  formulas  as  M  and  N  in  the  FEM-C 
model,  equations  (2-56)  and  (2-57).   However,  the  eigenfunc- 
tions  $ • (Z)  are  not  defined  at  the  same  levels  in  both  models. 

r\  vi" 

The  matrix  Q  is  defined  for  terms  of  the  form  -ttt'V, 
=  dZ 


i+1    i+1         T  d<})  . 
J  =  i-1  k=i-l   J  Zn 


The  matrix  P  is  defined  for  terms  of  the  form  ^tt'W, 

=  dZ 


2 

i+1  '  i+1  T  dcf>  . 

1=1  ^   Ui   /   dz1  \^idz  *  {2"73) 

j=i-l  k=i-2   D  Z     aL        K  1 


The  matrix  R  is  defined  for  terms  of  the  form  -^y'W, 


2 

i+1    i+1  T  dcj)  . 

R   =     I  I         T     j  ^   ^    dZ    .  (2-74) 

j=i-l  k=i-2   :  ZQ  aL        K  1 


The  staggered  eigenfunctions  present  three  general 
problems  for  evaluating  the  elements  of  the  five  matrices. 
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First,  for  an  N-layer  model,  portions  of  eigenf unctions 
(J>1  (Z)  and  <t>N+2(z)  are  defined  in  the  model  atmosphere  but 
the  physical  meaning  of  contributions  from  those  terms  is 
unclear. "  The  contributions  are  included  in  the  first  two 
rows  and  the  last  row  of  each  matrix.   Second,  only  portions 
of  eigenf  unctions  <J>2(Z)  and  <j>   ,  (Z)  are  defined  in  the  model 
atmosphere.   To  describe  the  incomplete  sides  of  both  eigen- 
functions  an  assumption  must  be  made  about  the  value  of  <J>2 
at  the  surface  and  $N+-,  at  the  top  of  the  atmosphere.   Last, 
the  equations  derived  for  the  general  elements  of  each  matrix 
are  more  complex  than  the  formulas  for  the  model  FEM-C.   The 
complexity  of  these  equations  makes  it  much  more  difficult 
to  evaluate  the  elements  in  the  first  two  and  last  two  rows 
of  each  matrix. 

Assumptions  are  made  and  procedures  are  developed  in 
an  attempt  to  resolve  these  problems.   In  this  model  the  mean 
state  variables,  u  and  T,  are  defined  only  at  the  N  staggered 
levels.   However,  u  and  T  values  defined  at  the  nodal  points 
of  <}>-j(Z)  and  <j>   „(Z)  are  important  in  the  Galerkin  formulation 
of  the  du/dZ  and  3T/8Z  terms.   In  my  thesis  experiments,  the 
values  of  u  and  T  are  defined  at  the  surface  and  top  of  the 
atmosphere;  they  are  not  defined  at  the  nodal  points  of  cj>,  (Z) 
and  cj)   _(Z)  becuase  these  points  are  outside  the  model  atmos- 
phere.  For  constant  shear  with  height,  u  and  T  are  defined 
at  the  boundaries  such  that  the  shear  in  the  two  half  layers 
at  the  boundaries  is  the  same  as  the  shear  in  the  other  layers 
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The  contributions  from  the  perturbation  quantities  defined 
at  the  nodal  points  of  <{>,(Z)  and  <J>N+2(Z)  are  not  included  in 
these  experiments.   To  be  consistent  with  the  other  N-layer 
models,  all  perturbation  quantities  are  defined  as  column 
vectors  of  length  N+l  even  though  all  quantities,  except  ver- 
tical velocity,  are  defined  at  only  N  levels.   The  values  of 
each  perturbation  at  its  lowest  staggered  level  are  stored 
as  the  second  element  in  its  column  vector;  zero  is  stored 
as  the  first  element,  which  corresponds  to  the  contribution 
from  (j)^(Z).   To  evaluate  the  staggered  eigenf unctions  defined 
in  the  layers  between  the  surface  and  Z~,  and  Z   ,  and  the 
top  of  the  atmosphere,  it  is  assumed  that  the  value  at  the 
boundaries  of  those  eigenf unctions  is  one-half.   Thus,  three- 
fourths  of  the  eigenfunctions  <{>2(Z)  and  ^.o^Z)  are  defined 
in  the  model  atmosphere. 

The  equations  for  the  general  elements  of  the  five 
matrices  are  evaluated  by  substituting  into  equations  (2-56)- 
(2-57)  and  (2-72)  -  (2-74)  the  formulas  for  <J>i  +  1(Z),  $i(Z), 
4>i_1(Z)  ,    ipi+1(Z),  ^i(Z),  ^i_1(Z),  and  ip.  2(Z)  terms  of  the 
local  coordinate  £  =    Z  -Z-.   The  equations  for  these  eigen- 
functions defined  for  the  levels  1,  2,  i,  and  N+l  are  listed 
in  Appendix  E.   The  equations  for  the  elements  of  the  first 
two  rows  and  the  last  row  of  the  matrices  P,  Q,  and  R  are 
derived  using  the  formulas  for  levels  1,  2,  and  N+l.   However, 
the  corresponding  rows  of  matrices  M  and  N  were  not  calcu- 
lated using  these  formulas  because  the  model  results  for  the 
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barotropic  experiment  were  not  constant  with  height  when  the 
elements  in  those  three  rows  of  both  M  and  N  were  calculated 
using  these  formulas.   The  reason  for  this  is  not  clear.   The 
equations  for  the  general  elements  of  M  and  N  are  simple 
enough  to  allow  the  terms  in  the  affected  three  rows  of  each 
matrix  to  be  deduced  by  reason.   The  equations  for  the  elements 
of  all  five  matrices  are  listed  in  Appendix  F. 

The  prognostic  matrix  equations  for  vorticity,  diver- 
gence and  temperature,  listed  in  Appendix  D,  are  simplified 
in  a  manner  similar  to  the  method  described  for  the-  FEM-C 
model.   The  final  form  of  the  forecast  equations  for  the 
vorticity  and  divergence  vectors  are  the  same  as  for  model 
FEM-C,  equations  (2-66) -  (2-69)  .   The  thermodynamic  vector 
equations  are 


Tl  _,,   =   T1BAR   ,  +  2At-[Ql-p-M   -N-T2 
— n+1      n-1         =—   H  =====  — 

+  (!■)  •M~1-Q-V1  -  M-1-R-W1]   ,  (2-75) 

K    ■     —    =====  —  n 


T2     ,  ,       =      T2BAR      ,     +    2At'[Q2+yM       -N-Tl 
— n+1  n-1  l^—      M   =========  — 

+    (4)M_1*Q-V2    -    M_1-R-W2]        .  (2-76) 

k  =====  —        =  —  n 


The  six  forecast  equations  are  solved  at  levels  two  to  N+1. 
The  surface  geopotential  and  the  diagnostic  variables  are 
calculated  using  the  corresponding  equations  in  model  FEM-A, 
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3.   FEM-B 

The  FEM-B  model  defines  vertical  velocity,  potential 
temperature,  mean  state  potential  temperature  and  diabatic 
heating  at  the  unstaggered  levels  in  terms  of  the  eigenfunc- 
tions  ty .  (Z) .   The  other  variables  are  defined  at  the  staggered 
levels  in  terms  of  the  eigenf unctions  <J>  •  (Z)  .   The  eigenf unc- 
tions are  the  same  as  defined  for  the  FEM-A  model,  shown  in 
Figure  3. 

The  finite  element  approximations  for  the  vorticity, 
divergence  and  thermodynamic  equations  are. derived  by  substi- 
tuting the  eigenfunction  expansion  for  each  dependent  variable 
into  equations  (2-35) -  (2-40)  .   The  vorticity  and  divergence 
equations  are  multiplied  by  <J>-  (Z)  and  integrated  with  respect 
to  Z  from  the  bottom  to  the  top  of  the  atmosphere.   The  resul- 
tant Galerkin  formulation  of  the  vorticity  and  divergence 
equations  are  the  same  as  those  derived  for  model  FEM-A.   The 
matrices  in  those  equations,  M,  N,  and  P,  are  the  same  as 
defined  for  FEM-A,  equations  (2-56) - (2-57)  and  (2-73).   The 
thermodynamic  equations  are  multiplied  by  ip- (Z)  because  poten- 
tial temperature  is  defined  at  the  unstaggered  levels.   As 
before,  the  equations  are  integrated  through  the  depth  of  the 
atmosphere.   The  resultant  equations  are  listed  in  Appendix 
G.   Four  additional  matrices  are  defined  for  the  two  thermo- 
dynamic equations.   The  mass  matrix  S  is 


i+1      ZT 

S   =     I  j         i>  .   (Z)ijj.  (Z)dZ  .  (2-77) 
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The  matrix  Q  is  defined  for  terms  of  the  form  grr'V, 


i+2    i+2        ZT 
Q   =     I  I         U.    /    <(>.  (Z)<J>.  (Z)i|;.  (Z)dZ  .  (2-78) 

j=i-l  k=i-l  J    Zn  J 


9T 
The  matrix  R  is  defined  for  terms  of  the  form  ~-s-*W, 


i+1    i+1        ZT  dip . 
R   =     J      IT./    gJ-  ^(Z)^i(Z)dZ  .  (2-79) 

j=i-l  k=i-l   J  Z- 


The  matrix  T  is  defined  for  terms  of  the  form  u«T, 


i+2    i+1        ZT 
T   =    I  7   u-    /    4>.  (Z)K  (Z)ij>.  (Z)dZ  .  (2-80) 

j=i-l  k=i-l   :  Z       J     K     ! 


As  discussed  in  the  FEM-A  model  description,  the 
staggered  finite  elements  present  problems  for  evaluating  the 
elements  of  the  matrices.   In  this  model  u  is  defined  at  the 
surface,  the  top  of  the  atmosphere,  and  at  the  N  staggered 
levels.   The  mean  state  temperature,  T,  is  defined  at  the 
unstaggered  levels  so  special  definitions  for  it  are  not  needed 
Also,  the  contributions  from  the  perturbation  quantities 
defined  at  the  nodal  points  of  (J),  (Z)  and  ty      5(Z)  are  not  in- 
cluded in  the  model  FEM-B  model  experiments.   The  staggered 
eigenf unctions,  cj>.(Z),  are  evaluated  at  the  boundaries  using 
the  assumptions  discussed  in  the  previous  section. 
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The  elements  of  matrices  S,  Q,  R,  and  T  are  evaluated 
by  substituting  formulas  for  <J>.  ~(Z),  <J>.  .(Z),  <J>  •  (Z)  ,  $.  ,  (Z)  , 
ij;.+1(Z)/  ^(Z)  and  4>j,(Z),  defined  in  terms  of  the  local 
coordinate  n  =  Z  -Z!f  into  equations  (2-77) -  (2-80) .   Formulas 
for  these  eigenf unctions  are  listed  in  Appendix  H.   The  equa- 
tions for  the  matrix  elements  are  listed  in  Appendix  I. 

The  forecast  matrix  equations  for  vorticity,  divergence 
and  temperature  are  simplified  in  a  manner  similar  to  the 
method  described  for  model  FEM-C.   The  final  form  of  the  vor- 
ticity and  divergence  vector  equations  are  the  same  as  for 
model  FEM-C,  equations  (2-66) -  (2-69)  .   The  thermodynamic  vec- 
tor equations  are 


Tl  J_-   =   T1BAR   ,  +  2At  [Ql  -  ji-S   -T«T2 

— n+1      n-1        ■=—    ==:  — 

+  (4)'S~1*Q-V1  -  S   -R-Wl]   ,  (2-81) 

k  =  —   =====  —  n 


T2  ,,   =   T2BAR   .  +  2At[Q2+yS  1T-Tl 
— n+1      n-1        —       ,  — 

. 
+  (^)'S_1-Q-V2  -  S"1-R-W2]   .  (2-82) 

K   ====  =====  n 


The  vorticity  and  divergence  equations  are  solved  at  levels 
two  to  N+1.   The  thermodynamic  equations  are  solved  at  all 
unstaggered  levels.   The  surface  geopotential  and  the  diagnos- 
tic variables  are  calculated  using  the  corresponding  equations 
in  model  FDM-B. 
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III.   EXPERIMENTS  AND  RESULTS 

Three  experiments  are  performed  with  each  model;  an 
initial  perturbation  in  the  meridional  flow,  flow  over  moun- 
tain topography  and  flow  with  a  diabatic  heat  source.   The 
first  two  experiments  are  performed  with  each  model  defined 
with  six  and  then  sixty  layers.   The  heating  experiment  is 
repeated  with  six,  twelve  and  sixty-layer  models.   The  analytic 
solution  of  each  experiment  has  not  been  derived.   For  each 
experiment,  the  60-layer  model  results  are  intercompared  to 
determine  if  the  models  are  converging  to  the  same  solution. 
The  standard  of  comparison  for  each  six-layer  model  is  its 
corresponding  60-layer  solution.   Temperature  and  divergence 
profiles  are  examined  in  each  experiment. 

Several  parameters  are  defined  identically  in  each  experi- 
ment.  The  vertical  coordinate,  Z,  is  defined  between  zero 
and  one  (1000-368  mb)  and  the  vertical  levels  are  equally 
spaced.   The  x-wavelength  is  4,000  kilometers.   The  time  step 
is  17.7  minutes.   The  Coriolis  parameter  is  defined  at  45 
degrees  latitude.   There  is  no  vertical  shear  in  the  u  field 
and  u  =  10.0  meters/second  (m/s) .   The  mean  state  potential 
temperature  increases  with  height  from  its  surface  value  of 
310.0  Kelvin  (K) . 

A.   ROSSBY  WAVE  EXPERIMENT 

Rossby  waves  are  generated  in  each  model  using  an  initial 
perturbation,  v*  =  5.0  m/s,  in  the  cosine  term  of  the 
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meridional  flow.   All  other  perturbations  are  initially  zero. 
There  is  no  diabatic  heat  source  and  no  mountain  topography. 
The  latitudinal  variation  of  the  Coriolis  parameter,  B,  is 
defined  at  45  degrees  latitude.   The  forecast  length  of  this 
experiment  is  9  6  hours. 

1.   Sixty-Layer  Models 

The  60-layer  FDM-A,  FDM-B,  FDM-C,  and  FEM-C  models 
converge  to  the  same  temperature  and  divergence  solutions 
(Figures  4-7) .   It  should  be  noted  that  the  phase  of  each 
variable  is  defined  between  zero  and  360  degrees.   There  is 
a  discontinuity  in  the  phase  profile  if  the  phase  passes 
through  zero  degrees.   The  height  at  which  the  temperature 
phase  discontinuity  in  model  FDM-A  occurs  differs  from  the 
other  three  models  because  temperature  is  defined  at  the 
staggered  levels  in  FDM-A.   The  four  models  represent  the 
same  physical  solution,  which  is  called  the  consensus  solution. 

The  FEM-B  temperature  amplitude  is  slightly  larger 
than  the  consensus  amplitude  (5%)  and  an  amplitude  oscillation 
is  present  in  the  lowest  three  layers  of  the  atmosphere  (Figure 
8) .   The  jagged  profile  may  be  caused  by  the  terms  in  the 
matrices  which  represent  contributions  from  the  eigenfunctions 
near  the  lower  boundary  of  the  model.   There  is  a  discontinuity 
between  the  temperature  phases  in  the  lowest  two  layers  of  the 
FEM-B  profile  (Figure  9).   Apart  from  this  feature,  the 
general  shape  of  the  phase  profile  is  similar  to  the  consensus, 
although  the  phase  is  five  degrees  less  and  passes  through 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (60  LAYERS) 


0.00  0.02  0.04  0.06  0.08  0.10 

TEMPERATURE  AMPLITUDE  (KELVIN) 


0.12 


0.14 


LEGEND 
□  =  FDM-A 
o  =  FDM-B 
a  =  FDM-C 
+  =  FEM-C 


Figure  4.   Sixty-layer  Rossby-wave  experiment  at  96 

hours.   Temperature  amplitude  profiles  are 
compared  for  models  FDM-A,  FDM-B,  FDM-C, 
and  FEM-C. 
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TEMPERATURE  PHASE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  5.   As  in  Fig.  4  but  for  temperature  phase. 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure 


6.   As  in  Fig.  4  but  for  divergence  amplitude 
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DIVERGENCE  PHASE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  7.   As  in  Fig.  4  but  for  divergence  phase 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  8.   Sixty-layer  Rossby  wave  experiment  at  9  6 

hours.   FEM-B  temperature  amplitude  profile 
is  compared  with  the  temperature  amplitude 
profile  of  FDM-C,  which  represents  the 
consensus  profile. 
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TEMPERATURE  PHASE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  9.   As  in  Fig.  8  but  for  temperature  phase. 
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zero  at  a  lower  level  than  the  consensus.   The  divergence 
amplitude  of  model  FEM-B  is  within  five  percent  of  the  con- 
sensus (Figure  10)  but  its  phase  profile  (Figure  11)  is  quite 
different  from  the  consensus. 

The  model  FEM-A  temperature  and  divergence  ampli- 
tudes vary  significantly  from  the  consensus.   The  temperature 
amplitude  profile  (Figure  12)  has  a  large  amplitude  oscilla- 
tion in  the  lowest  three  levels  above  the  surface,  possibly 
caused  by  the  choice  of  eigenfunction  representations  used 
to  evaluate  elements  of  the  first  two  rows  of  the  matrices. 
The  FEM-A  temperature  amplitude  is  larger  than  the  consensus 
throughout  the  atmosphere.   The  temperature  phase  passes 
through  zero  at  a  lower  level  and  is  five  degrees  less  than 
the  consensus  (Figure  13).   The  maximum  FEM-A  divergence 
amplitude  occurs  near  the  middle  of  the  atmosphere  and  its 
magnitude  is  ten  percent  less  than  the  maximum  of  the  consensus 
(Figure  14).   The  consensus  divergence  phase  is  nearly  con- 
stant with  height,  but  the  phase  profile  for  model  FEM-A 
steadily  decreases  with  height  (Figure  15) . 
2 .   Six-Layer  Models 

The  comparison  of  six  and  sixty-layer  profiles  for 
variables  defined  at  staggered  levels  may  be  initially  mis- 
leading.  The  first  staggered  level  in  a  six-layer  model  occurs 
at  Z  =  0.0833.   The  lowest  staggered  level  in  a  60-layer 
model  is  defined  at  Z  =  0.0083,  which  may  be  mistaken  for  the 
surface  in  the  graphs.   When  the  values  of  a  six-layer  model 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  10.   As  in  Fig.  8  but  for  divergence  amplitude 
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DIVERGENCE  PHASE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  11.   As  in  Fig.  8  but  for  divergence  phase 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  12.   As  in  Fig.  8  but  for  model  FEM-A  compared 
with  the  consensus,  FDM-C. 
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TEMPERATURE  PHASE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  13.   As  in  Fig.  12  but  for  temperature  pha 


se 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  14.   As  in  Fig.  12  but  for  divergence  amplitude 
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DIVERGENCE  PHASE  FOR  V  =  5.0  CASE  (60  LAYERS) 
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Figure  15.   As  in  Fig.  12  but  for  divergence  phase 
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coincide  with  a  60-layer  profile,  the  models  are  considered 
to  represent  the  same  physical  solution,  even  though  the 
six-layer  model  has  a  smaller  vertical  domain  for  staggered 
variables. 

The  temperature  amplitude  profiles  of  the  six-layer 
FDM-A,  FDM-B,  FDM-C  and  FEM-C  models  are  identical  with  their 
corresponding  60-layer  results  (Figures  16-19).   The  previously 
discussed  problems  in  the  lowest  layers  of  models  FEM-A  and 
FEM-B  are  quite  evident  in  their  six-layer  profiles.   The 
temperature  amplitude  of  FEM-B  is  virtually  identical  with 
its  60-layer  profile  above  the  first  two  layers  (Figure  20) . 
The  six-layer  FEM-A  profile  (Figure  21)  does  not  agree  with 
its  60-layer  profile  below  Z  =  0.50  (600  mb) .   The  large 
amplitude  oscillation  in  this  profile  destroys  the  integrity 
of  the  solution  in  a  significant  portion  of  the  atmosphere. 

The  six-layer  divergence  amplitude  profiles  for  the 
grid  C  models  are  identical  with  each  other  (Figures  22-23), 
and  the  divergence  profiles  for  the  grid  B  models  are  nearly 
the  same  (Figures  24-25) .   All  four  models  approximate  their 
corresponding  60-layer  solutions  with  similar  accuracy.   The 
six-layer  FDM-A  model  does  not  approximate  the  curvature  of 
the  consensus  60-layer  divergence  profile  as  well  as  the 
grid  B  and  C  models  (Figure  26) .   The  six-layer  FEM-A  model 
has  difficulty  approximating  the  location  and  magnitude  of 
the  mid-atmosphere  divergence  maximum  of  its  60-layer  solution 
(Figure  27) . 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (FDM-A) 
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Figure  16.   Six-layer  Rossby  wave  experiment  at  96 

hours.   Temperature  amplitude  profiles  are 
compared  for  the  six-layer  and  60-layer 
FDM-A  models. 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (FDM-B) 
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Figure  17.   As  in  Fig.  16  but  for  model  FDM-B. 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (FDM-C) 
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Figure  18.   As  in  Fig.  16  but  for  model  FDM-C. 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (FEM-C) 
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Figure  19.   As  in  Fig.  16  but  for  model  FEM-C. 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (FEM-B) 
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Figure  20.   As  in  Fig.  16  but  for  model  FEM-B. 
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TEMPERATURE  AMPLITUDE  FOR  V  =  5.0  CASE  (FEM-A) 
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Figure  21.   As  in  Fig.  16  but  for  model  FEM-A. 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (FDM-C) 
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Figure  22.   As  in  Fig.  16  but  divergence  amplitude  of 
FDM-C. 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (FEM-C) 
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Figure  23.   As  in  Fig.  22  but  for  model  FEM-C 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (FDM-B) 
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Figure  24.   As  in  Fig.  22  but  for  model  FDM-B. 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (FEM-B) 
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Figure  25.   As  in  Fig.  22  but  for  model  FEM-B 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (FDM-A) 
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Figure  26.   As  in  Fig.  22  but  for  model  FDM-A. 
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DIVERGENCE  AMPLITUDE  FOR  V  =  5.0  CASE  (FEM-A) 
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Figure  27.   As  in  Fig.  22  but  for  model  FEM-A. 
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B.   MOUNTAIN  TOPOGRAPHY  EXPERIMENT 

The  forced  vertical  velocity  term,  MTS,  in  the  surface 
geopotential  forecast  equation  (2-3)  is  non-zero  in  this 
experiment.   It  represents  the  contribution  to  surface  geo- 
potential from  air  flowing  over  mountain  topography  which 
varies  sinusoidally  in  the  x-direction  and  has  no  variation 
in  the  y-direction.   The  mountain  ridge-to-valley  height 
difference  is  1,500  meters.   To  reduce  the  trauma  for  the 
model,  the  mountains  are  gradually  "built"  to  their  full 
height  over  a  period  of  36  hours.   Thus,  the  forced  vertical 
velocity  increases  in  the  first  36  hours  of  the  forecast 
period  and  is  constant  for  the  remainder  of  the  96-hour  fore- 
cast period.   The  equations  used  to  define  the  forced  verti- 
cal velocity  are  included  in  Appendix  J.   Beta  and  all  initial 
perturbations  are  zero  in  this  experiment. 

1.   Sixty-Layer  Models 

The  FDM-A,  FDM-B,  FDM-C,  and  FEM-C  models  converge 
to  the  same  physical  solution.   The  amplitude  and  phase  pro- 
files of  temperature  and  divergence  are  depicted  in  Figures 
28-31. 

The  FEM-B  model  again  has  a  jagged  temperature  ampli- 
tude profile  in  the  lowest  two  layers  (Figure  32).   The 
temperature  amplitudes  in  the  remainder  of  the  atmosphere  are 
within  5%  of  the  consensus.   The  FEM-B  temperature  phase  has 
oscillations  in  the  bottom  three  layers  and  the  top  two  layers 
of  the  profile  (Figure  33) ,  but  the  rest  of  the  profile  is 
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Figure  28. 


Sixty-layer  mountain  topography  experiment 
at  96  hours.   Temperature  amplitude  profiles 
are  compared  for  models  FDM-A,  FDM-B,  FDM-C, 
and  FEM-C. 
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TEMPERATURE  PHASE  FOR  MOUNTAIN  CASE  (60  LAYERS) 
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Figure  29.   As  in  Fig.  2  8  but  for  temperature  phase 
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Figure  30.   As  in  Fig.  2  8  but  for  divergence  amplitude. 
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DIVERGENCE  PHASE  FOR  MOUNTAIN  CASE  (60  LAYERS) 
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Figure 


31.   As  in  Fig.  28  but  for  divergence  phase. 
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Figure  32.   Sixty-layer  mountain  topography  experiment 
at  96  hours.   FEM-B  temperature  amplitude 
profile  is  compared  with  the  temperature 
amplitude  profile  of  FDM-C,  which  represents 
the  consensus  profile. 
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Figure  33.   As  in  Fig.  32  but  for  temperature  phase. 
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within  one  degree  of  the  consensus  phase.   The  FEM-B 
divergence  amplitude  and  phase  are  identical  with  the 
consensus  values. 

An  oscillation  also  exists  in  the  lowest  three  layers 
of  the  temperature  amplitude  profile  for  model  FEM-A  (Figure 
34) .   The  amplitude  is  30%  higher  than  the  consensus  near  the 
surface, and  this  difference  decreases  with  height.   The  tem- 
perature phases  are  nearly  identical.   The  FEM-A  divergence 
amplitude  profile  is  similar  to  the  consensus,  but  the  amplitude 
does  not  decrease  as  fast  with  height  as  the  consensus 
(Figure  35) .   The  divergence  phase  profile  for  this  model  is 
the  same  as  the  consensus. 
2 .   Six-Layer  Models 

The  temperature  amplitude  profiles  of  the  six-layer 
models  FDM-A,  FDM-B,  FDM-C  and  FEM-C  are  identical  with  each 
other  and  also  with  the  consensus  solution  (Figures  36-39) . 
The  six-layer  staggered  finite  element  models  are  nearly 
identical  with  their  60-layer  solutions  above  the  first  two 
layers  (Figures  40-41) .   The  previously  identified  inadequate 
representations  of  the  temperature  amplitude  in  the  lowest 
layers  exist  in  this  experiment.   All  six  models  approximate 
their  60-layer  divergence  amplitude  profiles  quite  well 
(Figures  42-46) ,  however  the  FEM-A  model  (Figure  47)  does 
not  approximate  the  curvature  in  the  lower  portion  of  its 
60-layer  profile  as  well  as  the  other  models. 
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Figure  34.   As  in  Fig.  32  but  for  model  FEM-A  compared 
with  FDM-C. 
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Figure  35.   As  in  Fig.  34  but  for  divergence  amplitude 
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TEMPERATURE  AMPLITUDE  FOR  MOUNTAIN  CASE  (FDM-A) 


0.5  1  1.5  2  2.5 

TEMPERATURE  AMPLITUDE  (KELVIN) 


LEGEND 
c  =  FDM-A  60  LAYERS 
o  =  FDM-A  6  LAYERS 


Figure  36. 


Six-layer  mountain  topography  experiment  at 
96  hours.   Temperature  amplitude  profiles 
are  compared  for  the  six-layer  and  60-layer 
FDM-A  models. 
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Figure  37.   As  in  Fig.  36  but  for  model  FDM-B 
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Figure  38.   As  in  Fig.  36  but  for  model  FDM-C 
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Figure  39.   As  in  Fig.  36  but  for  model  FEM-C. 
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TEMPERATURE  AMPLITUDE  FOR  MOUNTAIN  CASE  (FEM-A) 
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Figure  40.   As  in  Fig.  36  but  for  model  FEM-A. 
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Figure  41.   As  in  Fig.  36  but  for  model  FEM-B 
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Figure  42.   As  in  Fig.  36  but  for  divergence  amplitude 
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Figure  43.   As  in  Fig.  36  but  divergence  amplitude 
for  model  FDM-B. 
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Figure  44.   As  in  Fig.  36  but  divergence  amplitude 
for  model  FDM-C. 
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DIVERGENCE  AMPLITUDE  FOR  MOUNTAIN  CASE  (FEM-B) 
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Figure  45.   As  in  Fig.  36  but  for  divergence  amplitude 
for  model  FEM-B. 
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Figure  46.   As  in  Fig.  36  but  for  divergence  amplitude 
for  model  FEM-C. 
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Figure  47.   As  in  Fig.  36  but  for  divergence  amplitude 
for  model  FEM-A. 
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C.   DIABATIC  HEATING  EXPERIMENT 

A  diabatic  heat  source  is  defined  in  the  layer  between 
Z  =  0.40  and  Z  =  0.60  (670-549  mb) .   The  rate  of  heating  is 
constant .in  time  and  varies  in  x  and  Z, 


2     (Z~  V 
Q(x,Z,t)       =      HEATING-cos     [-^= 4 — J  •  tt]  •  cos  (yx)     ,  (3-1) 

(ZU  _ZL) 


where  HEATING  is  5.0  K/day,  ZM  is  the  midpoint  of  the  heated 
layer,  and  Z_  and  Z„  are  the  lower  and  upper  boundaries, 
respectively,  of  the  heated  layer.   The  diabatic  heating 
vectors,  Ql  and  Q2_,  are  defined  in  the  initialization  subrou- 
tine and  stored  for  use  in  the  forecast  subroutine.   Beta  and 
all  initial  perturbations  are  zero.   The  forecast  length  is 
12  hours. 

1 .   Sixty-Layer  Models 

For  the  diabatic  heating  function  defined  in  equation 
(3-1),  the  maximum  heating  occurs  at  Z  =  0.50,  the  midpoint 
of  the  heated  layer.   The  models  defined  using  grids  B  and  C 
define  temperature  and  the  heating  functions  at  this  point. 
The  grid  A  models  do  not  have  temperature  and  diabatic  heating 
defined  at  this  point  so  the  maximum  rate  of  heating  in  these 
models  is  slightly  less  than  in  the  other  models,  and  the 
maximum  heating  occurs  throughout  one  layer  rather  than 
occurring  at  one  point.   The  heating  rate  at  each  level  is 
listed  in  Appendix  K  for  the  six,  twelve  and  sixty-layer 
staggered  and  unstaggered  models. 
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The  60-layer  profiles  for  the  six  models  are  quite 
similar,  but  the  differences  occur  because  the  models  are 
responding  to  different  forcing.   The  temperature  amplitude 
profiles -for  the  B  and  C  grids  come  to  a  sharp  point  at 
Z  =  0.50  and  the  grid  A  models  have  a  square-nosed  profile 
around  this  point  (Figures  48-50) .   The  height  in  the  atmos- 
phere  at;  which  the  model  experiences  the  onset  of  forcing  is 
slightly  different  between  the  grids.   This  explains  the  small- 
scale  oscillation  in  the  temperature  amplitude  profiles  near 
the  boundaries  of  the  heated  layer,  Z  =  0.40  and  Z  =  0.60. 
The  previously  identified  temperature  amplitude  oscillations 
in  the  lowest  layers  of  models  FEM-A  and  FEM-B  are  not  evident 
in  this  experiment.   The  temperature  phase  profiles  are  nearly 
identical,  except  that  the  grid  A  models  do  not  have  as  sharp 
of  a  spike  at  the  boundaries  of  the  heated  layer  as  the  other 
grids  (Figures  51-53) .   In  summary,  the  60-layer  temperature 
profiles  of  all  six  models  represent  the  same  physical  response 
to  the  diabatic  heating. 

Divergence  is  defined  at  the  midpoint  of  the  heated 
layer  only  in  the  grid  C  models.   Consequently,  the  shape  of 
the  divergence  amplitude  profile  for  the  grid  C  models  is 
different  than  the  other  four  models.   The  grid  C  profiles 
have  a  sharp  point  at  Z  =  0.50  (Figure  54) .   However,  the 
minimum  divergence  is  not  symmetric  around  Z  =  0.50  in  the 
grid  A  and  B  models  (Figures  55-57) .   The  divergence  ampli- 
tude is  identical  outside  the  heated  layer  for  all  models 
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Figure  48. 


Sixty-layer  diabatic  heating  experiment  at 
12  hours.   Temperature  amplitude  profiles 
are  compared  for  models  FDM-A,  FDM-B,  FDM-C, 
and  FEM-C. 
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Figure  49. 


Sixty-layer  diabatic  heating  experiment  at 
12  hours.   FEM-A  temperature  amplitude  profile 
is  compared  with  the  temperature  amplitude 
profile  of  FDM-C,  which  represents  the 
consensus  profile. 
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Figure  50.   As  in  Fig.  4  9  but  for  model  FEM-B  compared 
with  the  consensus,  FDM-C. 
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Figure  51.   As  in  Fig.  4  8  but  for  temperature  phase 
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Figure  52.   As  in  Fig.  49  but  for  temperature  phase 
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Figure  53.   As  in  Fig.  50  but  for  temperature  phase 
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Figure  54.   Sixty-layer  diabatic  heating  experiment  at 
12  hours.   Temperature  amplitude  models  are 
compared  for  the  two  grid  C  models. 
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Figure  55. 


Sixty-layer  diabatic  heating  experiment  at 
12  hours.   Temperature  amplitude  profiles 
are  compared  for  models  FDM-A  and  FDM-B. 
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Figure  56.   As  in  Fig.  55  but  models  FEM-B  and  FDM-C 
are  compared. 
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Figure  57.   As  in  Fig.  55  but  models  FEM-A  and  FDM-C 
are  compared. 
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except  FEM-A  (Figure  57) .   The  divergence  phase  profiles  are 
virtually  identical  for  models  FDM-A,  FDM-B,  FDM-C  and  FEM-C 
(Figure  58) .   The  divergence  phase  passes  through  zero  at  a 
different  level  for  model  FEM-B  than  the  consensus,  but  this 
difference  is  not  physically  significant  (Figure  59) .   The 
FEM-A  profile  differs  slightly  from  the  consensus  outside  the 
heated  }.ayer  (Figure  60)  .   In  conclusion,  the  six  models  are 
converging  to  the  same  60-layer  divergence  solution  inside 
the  heated  layer,  but  model  FEM-A  does  not  agree  with  the 
solution  of  the  other  five  models  outside  this  region. 
2 .   Six  and  Twelve-Layer  Models 

The  difference  between  grids  is  more  evident  in  this 
experiment  than  in  the  other  experiments.  Each  model  is  run 
with  both  six  and  twelve  layers  and  the  results  are  compared 
with  the  60-layer  consensus  of  the  six  models.  The  tempera- 
ture amplitude  consensus  profile  has  a  small-scale  decrease 
at  each  boundary  of  the  heated  layer  that  cannot  be  reproduced 
using  either  six-  or  twelve-layer  resolution. 

The  six-layer  grid  A  model  temperature  and  divergence 
fields  barely  respond  to  the  diabatic  heating.   The  grid  A 
models  have  identical  temperature  responses  for  both  six  and 
twelve  layers.   The  six-layer  perturbation,  constant  with 
height  in  the  heated  layer,  is  an  order  of  magnitude  smaller 
than  maximum  consensus  perturbation  (Figures  61-62).   The 
grid  A  models  have  a  stronger  response  to  the  heating  with 
twelve  layers  than  with  six  layers,  but  the  maximum  amplitude 
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Figure  58.   As  in  Fig.  48  but  for  divergence  phase 
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Figure  59.   As  in  Fig.  56  but  for  divergence  phase 
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Figure  60.   As  in  Fig.  57  but  for  divergence  phase 
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Figure  61.   Six-layer  diabatic  heating  experiment  at  12 
hours.   Temperature  amplitude  profiles  are 
compared  for  the  six-layer  and  60-layer 
FDM-A  models. 
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Figure  62.   As  in  Fig.  61  but  for  model  FEM-A, 
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is  only  60%  of  the  60-layer  consensus  response  (Figures 
63-64).   For  both  the  six-  and  twelve-layer  profiles,  the 
temperature  amplitude  at  the  two  staggered  nodal  points 
defined  closest  to  the  heated  layer  lie  directly  on  the  60- 
layer  profile.   This  may  indicate  that  for  much  less  than 
60-layer  resolution  the  amplitude  of  the  spike  will  not  equal 
the  maximum  amplitude  in  the  60-layer  experiment.   The 
divergence  fields  of  the  six-layer  grid  A  models  is  also  an 
order  of  magnitude  lower  than  their  corresponding  60-layer 
solutions  (Figures  65-66).   The  characteristic  sharp  divergence 
decrease  in  the  heated  layer  is  roughly  approximated  by  the 
12-layer  models  (Figures  67-68) .   The  minimum  divergence  is 
much  larger  than  in  the  consensus  profile  and  the  base  of  the 
spike  decrease  is  much  broader  than  in  the  60-layer  solutions. 
Outside  the  heated  layer,  the  shape  of  the  12-layer  profile 
is  similar  to  the  consensus,  but  the  amplitude  is  lower. 

The  six-layer  temperature  amplitude  response  of  the 
B  and  C  grid  models  are  virtually  identical  (Figures  69-72) . 
The  base  of  the  profile  spike  is  broader  in  the  six-layer 
models  than  in  the  consensus,  but  the  four  models  closely 
approximate  the  magnitude  and  width  of  the  tip  of  the  spike. 
The  12-layer  temperature  profiles  for  grids  B  and  C  are  also 
identical  (Figures  73-76) .   The  profiles  slightly  underesti- 
mate the  magnitude  of  the  perturbation  in  the  top  half  of  the 
spikes,  but  the  width  at  the  base  of  the  spikes  are  the  same 
as  in  the  60-layer  consensus.   The  small  scale  dip  at  the 


107 


TEMPERATURE  AMPLITUDE  FOR  HEATING  CASE  (FDM-A) 


-   'TP 


0.25  0.50  0.75  1  1.25  1.50  1.75  2 

TEMPERATURE  AMPLITUDE  (KELVIN) 


2.25 


LEGEND 
o  =  FDM-A  60  LAYERS 
c  =  FDM-A  12  LAYERS 


Figure  63. 


Twelve-layer  diabatic  heating  experiment  at 
12  hours.   Temperature  amplitude  profiles 
are  compared  for  the  12-layer  and  60-layer 
FDM-A  models. 
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Figure  64.   As  in  Fig.  63  but  for  model  FEM-A. 
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Figure  65.   As  in  Fig.  61  but  for  divergence  amplitude 
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Figure  66.   As  in  Fig.  62  but  for  divergence  amplitude 
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Figure  67..   As  in  Fig.  63  but  for  divergence  amplitude 
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Figure  68.   As  in  Fig.  64  but  for  divergence  amplitude 
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Figure  69.   As  in  Fig.  61  but  for  model  FDM-B 
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Figure  70.   As  in  Fig.  61  but  for  model  FEM-B 
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Figure  71.   As  in  Fig.  61  but  for  model  FDM-C. 
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Figure  72.   As  in  Fig.  61  but  for  model  FEM-C. 
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Figure  73.   As  in  Fig.  63  but  for  model  FDM-B 
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Figure  74.   As  in  Fig.  63  but  for  model  FEM-B. 
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Figure  75.   As  in  Fig.  6  3  but  for  model  FDM-C 
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Figure  76.   As  in  Fig.  63  but  for  model  FEM-C 
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boundaries  of  the  heated  layer  in  the  consensus  profile  is 
not  reproduced  by  the  12-layer  models,  otherwise  the  profiles 
agree  with  the  consensus  outside  the  heated  layer. 

The  diabatic  heating  causes  different  divergence 
amplitude  responses  for  the  B  and  C  grids.   The  grid  C  models 
have  similar  six-level  profiles  (Figures  77-78) .   The  minimum 
divergence,  which  occurs  at  the  center  of  the  heated  layer,  is 
nearly  the  same  in  the  six  and  sixty-layer  models.   There  are 
differences  between  the  two  grid  C  models  away  from  the  diver- 
gence minimum.   Outside  the  heated  layer  the  six-layer 
divergence  amplitude  of  the  FEM-C  model  provides  the  closer 
approximation  of  the  60-layer  consensus  profile.   Similarly, 
the  12-layer  FEM-C  profile  (Figure  79)  is  a  better  approxima- 
tion of  the  consensus  than  the  12-layer  FDM-C  model  (Figure 
80)  . 

For  both  six  and  twelve  layers,  the  grid  B  models  very 
poorly  approximate  the  60-layer  divergence  amplitude  within 
the  heated  layer.   The  12-layer  profiles  (Figures  81-82) 
are  identical  for  the  grid  B  models  and  are  closer  approxima- 
tions of  the  consensus  60-layer  profile  than  the  six-layer 
results.   The  12-layer  model  is  asymmetric  within  the  heated 
layer  and  does  not  have  the  dramatic  decrease  in  divergence 
that  exists  in  the  consensus.   The  twelve  and  sixty-layer 
profiles  are  similar  outside  the  heated  layer. 
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Figure  77.   As  in  Fig.  61  but  divergence  amplitude  for 
model  FDM-C. 
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Figure  78.   As  in  Fig.  61  but  divergence  amplitude  for 
model  FEM-C. 
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Figure  79.   As  in  Fig.  63  but  divergence  amplitude 
for  model  FEM-C. 
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Figure  80.   As  in  Fig.  63  but  divergence  amplitude 
for  model  FDM-C. 
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Figure  81.   As  in  Fig.  63  but  divergence  amplitude  for 
model  FDM-B. 
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DIVERGENCE  AMPLITUDE  FOR  HEATING  CASE  (FEM-B) 
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Figure  82.   As  in  Fig.  63  but  divergence  amplitude  for 
model  FEM-B. 
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IV.   CONCLUSIONS 

Although  the  analytic  solutions  of  each  experiment  have 
not  been  derived,  some  conclusions  about  the  characteristics 
of  each  model  can  be  drawn.   During  the  Rossby  wave  and  moun- 
tain topography  experiments  the  staggered  finite  element 
models,  FEM-A  and  FEM-B,  display  unusual  temperature  ampli- 
tude behavior  in  the  lowest  two  layers  of  the  model.   The 
oscillation  in  the  temperature  profiles  of  both  models  may 
be  generated  by  the  matrix  elements  which  represent  the 
contributions  from  the  eigenf unctions  near  the  surface.   It 
is  possible  that  better  representations  of  the  staggered 
eigenf unctions  near  the  lower  boundary  may  reduce  or  eliminate 
this  problem.   Jagged  temperature  profiles  were  not  observed 
in  the  diabatic  heating  experiment.   The  12-hour  forecast 
period  in  the  heating  experiment  may  be  too  short  for  the 
profile  discontinuities  to  grow  to  substantial  amplitudes. 

The  differences  between  the  grids  is  most  apparent  in  the 
diabatic  heating  experiment.   The  60-layer  temperature  pro- 
files of  all  six  models  represent  the  same  physical  solution. 
The  FEM-A  model  has  a  slightly  different  60-layer  divergence 
response  outside  the  heated  layer  than  the  other  five  models. 
The  grid  B  and  C  models  have  identical  temperature  amplitude 
profiles  for  both  six  and  twelve  layers.   However,  the  diver- 
gence amplitude  profiles  are  quite  different  and  the  grid  C 
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models  provide  the  better  approximation  for  both  resolutions. 
The  grid  A  models  barely  respond  to  the  diabatic  heating  with 
six  layers,  and  poorly  approximate  their  consensus  profiles 
with  twelve  layers.   There  are  indications  that  both  grid  A 
models  may  not  converge  to  their  60-layer  temperature  and 
divergence  solutions  with  much  less  than  60-layer  resolution. 

In  all  experiments,  the  grid  B  model  produces  virtually 
identical  results,  aside  from  the  small  amplitude  oscillation 
in  the  FEM-B  temperature  profiles.   Based  on  these  experiments, 
no  accuracy  is  gained  by  using  the  finite  element  approxima- 
tions with  grid  B.   In  two  of  the  three  experiments  the  grid 
C  results  are  identical.   The  FEM-C  model  provides  a  closer 
approximation  of  the  consensus  60-layer  divergence  profile 
than  the  FDM-C  model  in  the  diabatic  heating  experiment. 
Generally,  the  lower  resolution  FEM-A  model  does  not  approxi- 
mate its  60-layer  solution  as  well  as  the  FDM-A  model.   It  is 
not  known  if  the  different  approximation  characteristics  of 
the  FEM-A  model  are  due  solely  to  its  profile  problems  at  the 
bottom  of  the  atmosphere. 

The  importance  of  the  <J)-,(Z)  and  <j>   ~(Z)  eigenf unctions  in 
the  staggered  finite  element  models  is  not  clear.   They  are 
prominent  terms  in  the  mathematical  derivations,  but  their 
physical  significance  is  not  apparent. 
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APPENDIX  A 
FINITE  DIFFERENCE  APPROXIMATIONS 


1.   For  terms  of  the  form  -r=-  W: 

dZ 

a.   FDM-A  and  FDM-B ,  at  level  Z  =  Z.: 


,— y       ,      u.1?-u.  u.  -  u.  , 

dZ        2L  i  ^z.^,  -Z.;     i-l  4.  -Z.  /J 

l+l    i  i    l-l 


b.   FDM-C,  at  level  Z  =  Z!: 


-,—       n       u..-u.      u.  -  u.  , 

dZ  W      2   l  UZ!  .  -  Z!;    lZ!  -  Z!  /  J 

l+l    l       l    l-l 


3T 
2.   For  terms  of  the  form  -^=-  W: 

dZ 

a.   FDM-A,  at  level  Z  =  Z,: 


sm        i      T-^i  "  T-  T.  --T.  . 


b.   FDM-B  and  FDM-C,  at  level  Z  =  Z.: 

l 


—  T      —  T  T   —  T 

3Z  W      2   i  UZ!  .  -Z!   '     lZ!  -Z!  ,  H 

l  +  l    l         l    l-l 


3.   For  terms  of  the  form  -=-=-  V: 

dZ 

a.   FDM-A,  at  level  Z  =  Z.: 

1 


a77  i  u-,-.~U.       U.-U., 

dZ        2  v  UZ.  ,  -Z.  '     KZ.    -  Z.  ,  n 

l+l    1       1    l-l 
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b.   FDM-B,  at  level  Z  =  Z!: 


dZ  V      2^  i+l  +  vi;  KZ.^   -Z. 

l+l    i 


c.   FDM-C,  at  level  Z  =  Z!: 


,—        ,      U..-U.      u.-u., 

1+1    1      1   1-1 
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APPENDIX  B 
GALERKIN  FORM  OF  FEM-C  PROGNOSTIC  EQUATIONS 

1.   Vorticity  Equations  (2-35) -  (2-36)  : 


i+1   dA^1     ZT  i+1  ZT 

1        'Trr-        I         cj).  (Z)<|>.  (Z)dZ   =   -f    J"    D^(t)    /   <j>.  (Z)<j>.  (Z)dZ 
j-i-l     Z0     J  j-i-l       ZQ     3 

i+1    i+1   _   .         ZT 
-  y    I  I        u.A„(t)    /   (J).  (Z)<f>,  (Z)0.  (Z)dZ 

i=i-l  k=i-l  J  z    z„    :    K    x 


i+1  ZT 


r3 

j=i-l   X    Z 


S    £    V^(t)    /   (J>.  (Z)<f)i(Z)dZ  (B-l) 


0 


i+1   dA^     ZT  i+1  ZT 

I    ^    /    <f>.  (Z)<J>.  (Z)dZ   =   -f        D^(t)    /   (J).  (Z)«J).  (Z)dZ 
j-i-l      Z0  3=i-l        ZQ     3 

i+1    i+1  _  k        T 
+  y    J      I    u.A,  (t)    /   <j>.  (Z)«j>,  (Z)c|>.  (Z)dZ 
j=i-l  k=i-l   :       ZQ     J  1 

i+1    ,        ZT 


J 
2 
j=i-l   z     Z 


-6    I    V^(t)    /   <J>.  (Z)c{)i(Z)dZ  (B-2) 


0 
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a.   Divergence  Equations  (2-37) -  (2-38)  : 


i+1   dD?     ZT  i+1   •        ZT 


I        ^        I      <j>.(Z)<j>.  (Z)dZ  =  f    I     hUt)         j      <|>.(Z)<fr.  <Z)dZ 
o  3-i  1  z0 


i+1    i+1  _   .        ZT 

-y   I     I  u.D9(t)   /  <j>.  (Z)$,  (z)<j>.  (z)dz 

j  =  i-l  k=i-l  J  ^     Z      -1  x 


i+1   •        ZT 
-6    J  UJ(t)    /   <f>.  (Z)<j>.  (Z)dZ 

j-i-l     z0    ' 


2 

i+1    i+1  _  .        T  d<j>  . 

-y   I    J  u,w  (t)   /   _l  $  (Z)<f>.  (zuz 

j=i-l  k=i-l  3    z  zn    dZ   K    1 


„   i+1   •    '    ZT 
+  y     [   H^(t)    /   $.  CZ)  4> (Z)dZ  (B-3) 

j-i-l       2Q     3 


i+1   dD^     ZT  i+1   ■        ZT 

I      -fir-      j     cj).  (z)4>i(z)dz  =  f   £  A^(t)   /  <j>.(z)<f>i(z)dz 

j-i-l      ZQ  j-i-l       zQ 

i+1    i+1     .         ZT 

+   V   J      J   u.D^(t)    /   <J>.  (Z)<J>.  (Z)(J).  (Z)dZ 

j-i-l  k=i-l  ]  L         zn 


o 


i+1  ZT 


-  8   I  U?(t)    /  <)>.  (Z)<j>.  (Z)dz 
i=i-l      z„    J 


i    ■  "0 


z. 

r 

j-i-l  k=i-l         Z, 


i+1    i+1     .         T  d0 . 
+  V  m     I  l       u.W1(t)    /    ggl  d>k(Z)4,.(Z)dZ 


0 


,2   ^J 


+  y     £   H^(t)    /   cj).  (Z)<J>.  (Z)dZ  (B-4) 

j=i-i         z0      ^ 
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3.   Thermodynamic  Equations  (2-39 )- (2-40) 


i+1   dT^     ZT 
I        7TT-        I      <t>,  (Z)<J>  (Z)dZ 
j-i-i  at-  zn    ■« 


•  -,--,  zm 

l+l    l+l  _   ,  T 


=   -  y    I  I      u-T^(t)    /   cj>.  (Z)<j>.  (Z)«j>.  (Z)dZ 

j=i-l  k=i-l  J  Z     ZQ     :     K     X 

-    i+1    i+1  _  k        T  d<J>. 
j=i-l  k=i-l  J       ZQ 

2 
i+1    i+1  J  .         T  d(j)  . 

I  ,  f_J  ,Tjwi(t)   /   dz1  4>k<z)<,,i(z)dz 


j=i-l  k=i-l  J  *     Z 


0 


i+1    •         ZT 
+    £    Q^(t)    /   <{>.  (Z)((>i(Z)dZ  (B-5) 

j-i-l       z0 


i+1   dTp     ZT 
I    3^-    /   <j)i(Z)(j>i(Z)dZ 
j-i-1      zQ     3 


z 

i+1    i+1  _   .         T 
=   y    J      I      u  T  (t)    /   $  .  (Z)d>.  (Z)  <f>.  (Z)dZ 
j=£-l  k=i-l  -1  x     ZQ     -1 

2 
f    i+1    i+1    .         T  d(|)  . 

+  '!»    ?,    I  /jV2(t»  ,  /    n1  ♦]<<»>♦!<«■« 
3=1-1  k=i-l  J       Z_ 

2 
i+1    i+1     .         T  df 

".I      I   TjWjW    /    ggl  *k<Z)*.(Z)dZ 

j=i-l  k=i-l  J       Z0 

i+1   •         ZT 
+    J   Q^  (t)    /   (j).  (Z)(J>.  (Z)dZ  (B-6) 

j=i"l       ZQ     3 
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APPENDIX  C 
MATRIX  ELEMENTS  FOR  FEM-C 


1.   Notation:   A.   =   Z.  -Z.  ,,  but  A~  =  Z~  -Z! 

1       1    l-l        2     2    l 


A!   =   Z!  -Z!  , 
l       l    l-l 


u.  ,  T.   =   u,  T  at  i   vertical  level 
li      1 


2.   Matrix  definition:   M  (row,  column) 

The  i    row  of  a  matrix  contains  the  weighting  factors  of 
the  terms  that  affect  the  value  of  the  forecast  variable  at 
the  i  '  vertical  level.   The  weighting  factor  in  the  j 
column  of  the  i    row  is  the  influence  on  the  forecast  varia- 
ble at  the  i    level  from  the  variable  at  the  j    vertical 
level.   Level  i  =  1  is  the  surface,  level  i  =  N+l  is  the  top 
of  the  atmosphere. 


3.   To  calculate  the  terms  of  each  matrix,  substitute  in  the 
matrix  formula  for  <J>.  ,(Z),  <$>■(%),    tj>._,(Z)  using  functions 
defined  in  local  coordinates  with  respect  to  <£.(Z).   Define 
£  =  Z  -Z|. 

The  eigenfunctions  in  local  coordinates  are 


5  +  A! 

~KT  -Ai  i  «  i  ° 

<j>,  (5)     =        I  X  r  (C-l) 


1 


-^+Al+l 
Al+l 


0  <  £  <  A!  . 
—  ^  —   i+. 
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1*1 


1-1 


f  '  AU1 


£=0 


f--A, 


<*> 


i#i 


<P 


<P 


i-1 


♦1+1(«) 


Ai+i 


0     <     c, 


-    "i+1     ' 


(C-2 


*i-l(^ 


-  r 


-A!     <    £    <    0     . 

l    —        — 


(C-3) 


To  calculate  the  matrix  elements  corresponding  to  the  surface 
and  the  top  of  the  atmosphere,  use  the  above  formulas  for 
i  =  1  and  i  =  N+l,  respectively. 

4.   Mass  matrix,  M: 

a.   General  formulas,  i  =  2,...,N: 


M(i,i-1) 


6  l 


(C-4) 


M(i,i)   =   j(Ai+1H-A.) 


:c-5 


M(i,i+1) 


b^i  +  l 


(C-6) 


b.   At  level  i  =  1 


M(l,l 


=   K 


:c-7) 
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M(l,2)   =   ^A2  (C-8) 


c.   At  level  i  =  N+l : 

M(N+1,N)   =   |an+1  (C-9) 


M(N+1,N+1)   =   jAN+1  (C-10) 


5.   Matrix  N: 

a.   General  formulas,  i  =  2,...,N 

N(i,i-1)      =     ^-Ai-  (u±  +ui_1)  (C-ll) 


N(i,i)      =     K-(Ai+l+Ai}    +  ir-(ui+l*Ai+l 


+   ui_1-Ai)  (C-12) 


N(ifi+1)       =      ijAi+1-  (ui+1  +ui)  (C-13) 


b.       At    level    i    =    1: 

N(l,l)       =      |  Uj^-Aj    +   yz  *VA2  (C-14) 


N(l,2)       =      X7*A2'  (u2  +ux)  (C-15) 


c.       At    level    i    =    N+l: 


N(i,i-1)       =      ^■•Ai(ui  +ui-;L)  (C-16) 


N(i,i)       =      j'Ai-ui    +   XY'Ai-ui_1  (C-17) 
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6.   Matrix  P: 

a.   General  Formulas,  i  =  2,...,N 

P(i,i-1)   =  |(ui-ui_1)  (C-18) 


P(i,i)   =  j(ui+1-ui_1)  (C-19) 


P(i,i+D   =   ^ui+1-ui)  (C-20) 


b.   At  level  i  =  1: 

P(l,l)   =   j(u2  -ux)  (C-21) 


P(l,2)   =   |(u2  -u±)  (C-22) 


c.   At  level  i  =  N+l: 

P(N+1,N)   =   |(uN+1-uN)  (C-23) 


P(N+1,N+1)   =   j(uN+1-uN)'  (C-24) 


7.   Matrix  R: 

a.   General  formulas,  i  =  2,...,N 

R(i,i-1)   =  ^(Tj_   -Ti_1)  (C-25) 


R(i,i)   =   |(Ti+1  -T±-1)  (C-26) 


R(i,i+1)   =   |(Ti+1  -Ti)  (C-27) 
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b.   At  level  i  =  1: 

R(l,l)   =   |(T2  -Tj)  (C-28) 


R(l,2)   =   |(T2  -Tx)  (C-29) 


c.   At  level  i  =  N+l: 

R(N+1,N)   =   ^(TN+]_  -TN)  (C-30) 


R (N+l, N+l)   =   |(TN+1  -TN)  (C-31) 
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APPENDIX  D 
GALERKIN  FORM  OF  FEM-A  PROGNOSTIC  EQUATIONS 

1.  The  vorticity  equations,  (2-35) -  (2-36)  ,  have  the  same 

form  as  the  vorticity  equations  in  FEM-C.   While  the  FEM-A 

equations  look  the  same  as  equations  (B-l)-(B-2)  in  Appendix 

B,  the  equations  for  <$>.  ,    cj>  .  and  cj>.  are  different  because  the 

1    J       K 

0(Z)  eigenfunctions  are  defined  for  the  staggered  levels  in 
model  FEM-A. 

2.  Divergence  Equations  (2-37) - (2-38) : 

i+1   dD^     ZT  i+1   .        ZT 

J        -^F-        j      4>i  (Z)«J).  (Z)dZ-=  f    I      A?(t)    /   (J).  (Z)<f>.  (Z)dZ 
j-i-l      2Q     3  j-£-l       ZQ     3 

•  -,  z 

1+1    l+l  _   k         T 
-  W    I      £   u.D  (t)    /   <J>.  (Z)(J),  (Z)<j>.  (Z)dZ 
j=i-l  k=i-l  ]  z  ZQ    3    K    1 

i+1   .       ZT 


-3   I   U^  (t)    /   <f>.  (Z)<|>.  (Z)dZ 
i=i-l      Z„    : 


i+1    i+1    ,         T  d(J)  . 
"U    [     J   uJ(t)    /   _1  4,  (z)(0.(Z)dZ 
j=i-l  k=i-l  3    Z  Z      dZ    K     i 


2   i+1   i        Zt 
+  y     J   H^1  (t)    /   (J,.  (Z)0,  (Z)dZ  (D-l) 

j-i-l     z0    : 
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i+1   dD^     ZT 


i+1 


I       dt~  I     ♦j(Z)*1(Z)dz  =  f  I     A^(t)   /  <|>.  (Z)<D±(z)dz 


j=i-l      Z 


j-i-l 


i+1    i+1 


-  _k 


+  U    I      I   u,D1(t)    /   (f)i(Z)(J)k(Z)<f)i(Z) 

-j=i-l  k=i-l  J       Zrt     J 


dZ 


i+1   .         T 


-  B    J  UJ(t)    /   <J).  (Z)<j)i(Z)dZ 


j-i-l 


i+1    i+1  _  .         T  d<|>  . 
j=i-l  k=i-l  J      Z_ 


i+1 


U     £   H:J(t)    /   ((>.  (Z)«J)i(Z)dZ 


j-i-l 


(D-2) 


3.   Thermodynamic  Equations  (2-39) -  (2-40)  : 


i+1   dT^     ZT 

I     at^     /    $i(z)$i(z)dz 


j-i-l     Z, 


i+1    i+1 


=  "  V    I 


-   _k 


J      J   u.T9(t)    /   <j>.  (Z)<MZ)$.  (Z)dZ 
j=i-l  k=i-l  :  ■     Zn    :         x 


-    i+1    i+1  _  .         T  d(J>  . 
+  (p)        I  ^v£(t)   /   ^  <j>k(z)<j)i(z)dz 

zo 


j=i-i  k=i-i  : 


i+1    i+1 


-  ,Jc 


'T  d<j>. 


J      I   Tjw1(t)    /    g-1  ^k(Z)4)i(Z)dZ 


j=i-l  k=i-l 


i  +  1   •         T 


+      nQl(t)    /   *j(Z)(j)i(Z)dZ 


j-i-l 


(D-3) 


142 


i  z 

i+1      dT^  T 

I        7GT-        !      <j)i(Z)4>i(Z)dZ 


j=i-l   ~w 


dt      „    '       TJ 
0 


i+1        i+1 


-  „k 


M-       I  I      u,T    (t)  /       (J).  (Z)(|>k(Z)((,i(Z) 

i=i-l   k=i-l    J  Zrt  J 


dZ 


i+1        i+1 


'T   d<j>  . 


+   (f)    .   I  l     u^Ct)       /       ^l^(z)^(z)dz 

1  =  1-1    k=i-l    J  Z„ 


i+1         i+1 


-  „k 


T   d<(>. 


I  I      Wit)         /        _1  ^k(Z)(|,i(Z)dZ 


j    2VW'     7    J  dZ 

j=i-l    k=i-l    J  Z 


i  +  1       •  T 


+        I      QUt)         j      <$>■  (Z)<b.  (Z)dZ 
i=i-l  Zrt  J 


(D-4) 
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APPENDIX  E 


LOCAL  COORDINATES  FOR  STAGGERED  (^(Z) 


1.   To  evaluate  the  integrals  in  the  finite  element  approxima- 
tions of  the  vorticity  and  divergence  equations  in  models 

FEM-A  and  FEM-B,  express  <$>•,-,'    <J>  •  i    $•_■]'  ^  ■ +±'    ^-'  ^-_i 

and  \b .    -,  in  terms  of  the  local  coordinate  E  =  Z  -Z-.   The 
ri-2  ^        1 

general  form  of  the  local  coordinate  system  is  shown  below. 


i+1 


A'. 


z 


*   2 


£  =  0 


5    2 


1-1 


f.-A, 


i-1 


i-1 


i-2 


2.   In  terms  of  local  coordinates  for  (j).  (£)  *  the  equations 

for  tt  .  .,  ,  6.,  <t>  ■    t/  V-  -i/  ^-*  ty-     i  and  'J; .  ~  are  listed  below, 
l+l    l    l-l    l+l    i    i-l       i-2 

The  notation  used  is  defined  in  Appendix  C. 
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5+A± 


<Pi(Z)      = 


-A.  <  C  <  0 
A .  l  —   — 

l 


.  1+x   0  <  C  <  A.^. 

A.+1       -    -   1+1 


♦i+i<5)    ■  r^T         °i5i  a1+1 

1+1 


W5)    =    if  -a±  1  C  1  0 

5-4*!  -A!       A! 


I 

A, 


-2"  <  5  1  2 


*i(5)   =        1,, 


-C  +tA'  +  A'  A! 

2  x i±i  _±  <  f    <  A 

Ai+1  2  -  «  -  Ai+1 

5-yA!  A! 

W^>  =  -TTT  "T^l  Ai+i 

i+i 

£  +tA-  +  A!  ,'  -A! 

2  l    l-l  ,   „  ,.  _    l 


-A.  <  C  < 


A!  ,  l  -  s  -  2 

l-l 


*i-i(5)   M   ..l.. 


^+2A1  "Ai    t    Ai 

AT"        -2-15  1  — 

1 


-5  -^A!  -A! 


1-1 


3.   With  respect  to  the  function  <j>.  {£,)  ,    the  layer  of  interest 
is  Zq  <_   Z  <_   Z2  or  yA'  1  C  1  A2 '   Eigenf unctions  $.,(£)  , 
<f>2(0/  <|>-,(£)  and  ip2(?)  are  defined  in  this  layer. 
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•0   i  i  i  i  i  r 


£  =  A  2 


0!  (O     = 


<MS>   = 


-C  +  A^ 


A  • 

"2 


<hU)     = 


*2(5) 


-?+|a- 

A2 
C-|a' 

-2 

4.   With  respect  to  the  function  $9(£)*  the  layer  of  interest 

-A  ' 
2 
is  ZQ  <_  Z  _<  Z,  or  -= —  <_•  £,    <_  A^.   Eigenf unctions  $,,  (J)-,  $_, 

ip,  ,  \p~    and  ^  are  defined  in  this  layer. 


f- A. 


i  i  i  i  i  i 


5  "  2 


£  =  0 


4    2 


*2 
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♦l<e>  =  if 


*2(C)   = 


„  _L 


*3(5) 


- 

|a*  <  5  1  0 

C+A' 
A2 

- 

A2 

-T<  5  <  0 

-5+A3 
A3 

0 

1  C  1  A3 

0 

1  ?  1  A3 

*XU)  = 


^+IA2 


*2^)   = 


A2 

/  «  4A2 
)     A2 

"A2       A2 
2   <  5  £  22 

)  "«  4A2 

v       a3 

+  A3 

A2 

-f  1  5  1  A3 

5-|a- 
A3 

A2 

*3(U   = 


5.   With  respect  to  the  function  <t>N,,  (O  /  the  layer  of 

interest  is  Z,T  <  Z  <  Z '  .  or  -A.T1,  <  E,    <    ~-A '  ,  ,  .   Eigenf  unctions 
N  —    —   N+l       N+l  —    —  2  N+l 

^N+2'  ^N+l'  *N'  ^N+l'  ^N  and  ^N-l  are  defined  in  this  layer. 
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■?'    lilt 


N+1 


ZN      — 


q  '    2 

A    m         — N*1 


f«-AN+1 


N-1 


«N+2(?) 


"N+l 


0     <    ^    < 


^+: 


-   •"   -        2 


*N+1(?> 


-£  +  A_' 


"N+l 


5  +^ 


N+l 


'N+l 


0   1   ?  i  2*8+1 


"Vi  i  «  i  ° 


*„(€> 


'N+l 


"AN+1    i    «   i   ° 


W(^        = 


--   +iA' 
2  "N+l 

~N+1 


1*. 

•   2^N+1 


0     <     C  T^ 


*N(«J 


s>       2^N+1 
.  i 

"N  +  l 

£  +— A'  +  A' 

^       2^N+1       "N 

/>  > 
N 


°^i   2-N+l 


■W    i    5    ±   ° 


»K-1(« 


2-N+l 

i     I 

JN 


- 1  • 


-A 


<     £     < 


'N+l 


N+l    -    s    -         2 
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APPENDIX  F 
MATRIX  ELEMENTS  FOR  FEM-A 

1.   Mass  matrix,  M,  using  the  notation  defined  in  Appendix  C: 
a.   General  formulas,  i  =  2,...,N: 

'   M(i,i-1)   =  |a±  (F-l) 


M(i,i)   =   j(Ai+1+Ai)  (F-2) 


M(i,i+1)   =   ^Ai+1  (F-3) 


b.   At  level  i  =  1: 

M(l,l)   =   jA2  (F-4) 


M(l,2)   =   |a2  (F-5) 


c.   At  level  i  =  N+l 

M(N+1,N)   =   ^-AN+1  (F-6) 


M(N+1,N+1)   =   i(AN+1  +  |A'+1)  (F-7) 


2.   Matrix  N: 

a.   General  formulas,  i  =  2,...,N 

N(i,i-1)   =  Yj(n±  +u±_1)  '&•_  (F-8) 
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N(i,i)       =      ^(U^-A.^+U.-A.) 


+  K-(Ai+i+Ai)  (F-9) 


N(i,i+1)      =      2T(ui+l  +ui}  '  Ai+1  (F-10) 


c.       At    level    i    =   N+l 


N 


N(i,i)       =      2T(VkTAn+1>     +    A(UN+l*Vl) 


3.   Matrix  P: 

a.   General  formulas,  i  =  3,...,N 


A. 
P(i,i-2)       =       (u.   -u 


i         i-l'     L6-A!     , 
l-l 

(A!_)3 


48-  (Ai)2-  (A|_x 


b.      At    level    i   =    1: 

N(l,l)       =      3^-(^2*A2)     +   ?(^1*A2)  (F-ll) 


N(l,2)       =      3^-(u2   +ux)  «A2  (F-12) 


(i,i-l)       =      i^N+l+V'Vn  (F'13) 


+  ivi^w+k+i1  (F-l4) 


(F-15) 
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Ai  (AP 

P(i,i-1)   =   {(u.  ,  -u.)  •  [g— r^ 2 

1+1    l    8  Ai+1    4g. (Ai+i)Z 


A! 

l 


Ai 


i-l 


4-AU 


8-  A. 


(A|) 


(A!) 


(Aj_)3 


48-  (Aj_) 


8*VAi-l 


48-  (Ai)  (A[_1 


■]  } 


(F-16) 


P(i,i)   = 


A! 

t(ui"ui-i),[n: 


(A!_) 


48-  (Aj_) 


3*Ai 

+    (Ui  +  1-Ui)-[gT^ 


5' (A|)2 


i+1    48*(Ai+1) 


+  (ui+l"ui: 


i     a: 
[- X 


+ 


A! 

l 


l  +  l 


2'Ai+l  '2'Ai+l   6*Ai+l 


(A!_) 


8'(Ai+1)(A!+1) 


(A1+1> 
B'(Ai+1) 


(a:; 


48.(Ai+1)  (A!+1) 


(F-17) 


P(i,i  +  D   = 


A! 

l 


A!) 


i+1 


4'Ai+l 


•Ai+1-Ai+1 


(A! 


48-(Ai+1)  (A!+1) 


b.   At  level  i  =  1: 


(F-18) 


P(l,l)   =   g(u2  -  u1 


(F-19 


P(l,2) 


4¥(u2  -Ul} 


(F-20) 
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c.   At  level  i  =  2: 


13  .-   - '.    .-   -  .  ,  A2      <A2»2 


P(2,l)   =   w-(u2  -Ul)  +  (u3  -u^-Ey—  -  ]        (F-21) 


'3    48-  (A3) 


5   -   -      -   -     1   A3    A2     A2 
P(2,2)   =  ^-(ua-u^  +  (u3-u2).[I-^  +  4fT-gT-A- 

(A^)2       (A2)2      (A2)3 


+  2 = +  =-5 ]  (F-22) 

48(A3)     8-A3-A3    48(A3)  (A') 

A3      A2      (A2)2        (A2)3 


p(2,3)   =   (u3  -u2)-[^  -  jTaT  +  8- A, -AA  "  ~  ,2/<t  J 

(F-23) 


,3     h    a3     o  a3  za3     48(A3)"(A^) 


c.   At  level  i  =  N+l : 

AN+1        (AN+1>3 


P(N+1,N-1)   =   ("Mrt-Vt^Tg-  ~  ,,    ,2,..  '  (F"24) 


N    48-<AN+l)   <AN> 


P(N+1,N)   =   t^(uN+2  -u    )  +  (UN+1-UN).[|  -  g^±| 

N 

,  AN+1     AN+1   ,     (AN+1»2 


AN    8-AN+l    48MAN+1)2 

(A',)2  (A'  ,)3 

N+l      +  N+l  ]}  (p_25) 


8(AN+1)2(AN)    48(AN+1)2(An' 


A ' 

P(N+1,N+1)   =   {^(Vo-Vi)  +  (U.,^  -  l\J  •  [ 


48  v  N+2    N+l     v  N+l    N   l8-A^, 

N+l 

N+1    M  (F-26) 


48*(AN+1)2 
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4.   Matrix  Q: 

a.   General  formulas,  i  =  2,...,N 

0(1,1-1)   =   |(u±  -ui_1)  (F-27) 


Q(i,i)   =   j(^i+1-ui_1)  (F-28) 


Q(i,i+1)   =   ^(ui+1-ui)  (F-29) 


b.   At  level  i  =  1: 

Q(l,l)   =   -j(u2  -u^  (F-30) 


Q(l,2)   =   |(u2  -u1)  (F-31) 


c.   At  level  i  =  N+l: 

Q(N+1,N)   =   |(^j+1-V  (F-32) 


Q(N+1,N+1)   =   |(UN+2"UN)  (F-33) 


5.   Matrix  R: 

The  formulas  for  the  elements  of  this  matrix  are  the  same 
as  those  for  matrix  P,  with  T  replacing  u  in  equations 
(F-15)-(F-26) . 
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APPENDIX  G 
GALERKIN  FORM  OF  FEM-B  PROGNOSTIC  EQUATIONS 

1.  The  vorticity  equations  have  the  same  form  as  the  vor- 
ticity  equations  for  model  FEM-C,  equations  (B-l)-(B-2). 
However,  the  equations  for  $-,,  <{>  •  ,  and  <Jk  are  different 
because  4>(Z)  is  defined  for  the  staggered  levels  in  FEM-B. 

2.  The  divergence  equations  are  the  same  as  the  divergence 
equations  in  model  FEM-A,  equations  (D-l)-(D-2). 

3.  Thermodynamic  Equations  (2-39) -  (2-40 )  : 


i  +  1   dT^1     ZT 

I     g^     /    ^.(z)^  (z)dz 

j-i-i     zQ    3 

i+2    i+1     ,         ZT 
-  -  y    I  I      u,T  (t)    /   (f»,(Z)^  (Z)<j,  (Z)dZ 

j=i-l  k=i-l  J      ZQ    J         x 

f    i+2    i+2  ZT  d<£  . 

+  <§>  .  I   ,  .  I  UjVjCt)    /   ^  *k(zH.  (Z)d2 
j=i-l  k=i-l  J       Zn 

i+1    i+1     .        ZT  di|/. 

"  .  ?  ,  v  J  /j"!'"'  ,  /    dz^k'2'*!12"32 
J  =  i-1  k=i-l  J       Z 

i+1   •         ZT 
+    I       0?(t)    /   ip-(Z)^ .  (Z)dZ  (G-l) 
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i+1      dT^  ZT 


l       or-       I     ^i(z)^i(z)dz 


j=i-l   ~to      Z 


dt      „    J       T3 
0 


i+2         i+1   _     .  T 

=      U         I  I      u.T,  (t)  /       4>  .  (Z)  ipfc  (Z)  ip.  (Z)dZ 

j=i-l  k=i-l   J    x  z  J  1 


i+2        i+2   _     k 


T  d<J> . 


]=i-l    k=i-l    J  Z- 


i+1        i+1  k  T   d\p 

l      J,  /^(t)       /      ^i*k(z)*i(z)dz 

1=1-1    k=i-l    J  Z_ 


j=i-l    "  Z 


J         dZ 

i/j  . 


i+1       •  ~T 

I     Q2{t)        I       ipj(Z)^i(z)dz 


(G-2) 
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APPENDIX  H 
"  LOCAL  COORDINATES  FOR  UNSTAGGERED  \j> .   (Z) 

1.   To  evaluate  the  integrals  in  finite  element  approximations 

of  the  thermodynamic  equations  in  model  FEM-B,  express  <p .    ~ , 

<b .    ,  ,  <b  .  ,    d> .  ,/  to .  ,  -  ,  \b .    and  \i> .    ,  in  terms  of  the  local 
vi+l   yi'  yi-l   ri+l   ri      ri-l 

coordiante  n  =  Z  - Z ! .   The  general  form  of  the  local 
coordinate  system  is  shown  below. 


z.'  1  =  b.^ 


.£.♦, 


z,+, , 

Z.'  >?  =  0 

.; —  ...f 

z; , >?  =  -  a'. 

1-1  i 


2.   In  terms  of  the  local  coordinates  for  V-(n),  the  equations 

for  <yi+2'    *i+i'  *-i'  *i-l'  ^i+1'  vi  and  ^i-1  are  llsted  below- 
The  notation  used  is  defined  in  Appendix  C. 
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<fr±(Tl)         = 


^  +  Ai+ki 


-*  4Ai+i 


1+1 


A!    <   n    <  -  tA- 
i    —        —        2    i 


1.,  1. 

2    i   —        —  2    l+l 


cj>i_1(n)     = 


-n  -|a! 


♦i+1<n>      = 


c})i+2(n)      = 


A. 

i 

( 

n4A! 

) 

Ai+1 

I 

-n  +|a:  + 

ai+i 

V 

4i+i 

n 

-k+i 

i+l 


A!    <  n   <  -|a! 


-_A.  <  n  <  _A.+1 


jA!+1  i  n  i  a:+1 


|a!+1  <  n  i  A!+1 


ipi(n)       = 


n  +  A!_ 

aT~ 

i 


-a:  5  n  i  o 


o   <   n   <   A!+. 


*i-i(n)    -    ir 


-a!    <   n    <_  0 


*i+l(n)       =      A! 


i+l 


o  1  n  £  A!+1 


3.   With  respect  to  the  function  ijj,  (n)  ,  the  layer  of  interest 
is  ZQ  <_  Z  <_   Z'  or  0  <_   n  <_  A'  .   Eigenf unctions  $,  ,  <|>2,  <j>_, 
ip   and  (|;„  are  defined  in  this  layer. 
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1  I  I  I  I  I  I 


1  -  0 


$2  (n 


-n  +  2a' 


o  r)   <--a: 

_  2  2 


*2(n) 


n  ^jA2 


n   2"2 
"2 


0  <  n  <  |a2 


|a2  <  n  <  A2 


<Mn)   = 


2  "2 


^A'  <  r,  £  A' 


^(n)   = 


-^  +"2 
A2 


0  £  i  <_    A2 


*2(n) 


n 
/>  i 
"2 


0  1  r;  I  A2 


4.   With  respect  to  the  function  ijj-lr:),  the  layer  of  interest 
is  Zq    Z    Z-,.   Eigenf unctions  $,  ,  i-,  »  <j>_,  <£  .  ,  ^,  ,  i|/2  and 
ij;..  are  defined  in  this  layer.   Formulas  for  <j>  -  ,  i^,  ,  \\j0    and 
^3  are  obtained  from  the  general  formulas  with  i  =  2. 
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Z3 


Zn       I  I  I  I  I  I  I  I  I  I 
0 


"I  =  A3 

A3 

1  =  0 

1  =  -  A, 


0X  (n) 


A2 


-A'  <  n  <  -  jA' 


<P2  (n) 


n  4A2 

A2 

-n  +  |a' 


-^2  1  n  1  "  IA2 


4A2  1  n  1  lA3 


43(n 


n  +  -A1 


-r,  +^2  +^3 


2L2  -  n  -  2^3 


-A '  <  n  <  A ' 
2~3  -  n  -   3 


5.   With  respect  to  the  function  i!K+-,   (n)  ,  the  layer  of  interest 

is  Z'  <  Z  <  Z',1  or  -A'  ,  <  n  <  0.   Eigenf  unctions  <b„ ,  „ , 
N  —    —   N+        N+l  —    —        r  N+2 

$XT,i  »  <t>,,/  U;.T ,  ..  and  U;.T  are  defined  in  this  layer. 

N+l   ^N    N+l       N  J 
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,/     ■■■'■■'■■' 

N+1 


1  -  0 


Z        __ 7,  =  -  An+1 

N+1  2 


»?  =  -  A', 


N  +  1 


N+2 


N+  1 


N  +  1 


♦„("> 


n  -  —A* 
1       2"N+1 


A 


N  +  1 


-"N+1    -    H   -         2LN+1 


$N+l(n) 


"n   4AN+1 

A  ' 
N+1 

n+AN+l+IAN 
AN+1 


TAN+1   ^1° 


•  a  '        <   n   <   -  —  A ' 

"N+1    -         -         2"N+1 


*N+2(n) 


n  +  2an+i 

""N+1 


— A  *  <     r,     <     0 

2^N+1    -         -    U 


Vn) 


"N+1 


'n+i  1  n  £  o 


*H+l(n) 


n   +  A  '      , 
N+1 
■  i 

"N  +  1 


A'        ,        <       r       <       0 

n+1  -  n  - 
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APPENDIX  I 
MATRIX  ELEMENTS  FOR  FEM-B 

1.  Matrices  M,  N,  and  P  are  the  same  as  those  defined  in 
model  FEM-A.   The  equations  for  their  elements  are  listed  in 
Appendix  F,  equations  (F-l) - (F-26) . 

2.  Mass  matrix  S  is  the  same  as  matrix  M  in  model  FEM-C. 
The  equations  for  the  elements  of  S  are  the  same  as  equa- 
tions (C-4)-(C-10)  in  Appendix  C. 

3.  Matrix  R  is  the  same  as  matrix  R  in  the  model  FEM-C. 
Equations  (C-25) - (C-31)  in  Appendix  C  define  the  elements  of 
R. 

4.  Matrix  Q,  using  the  notation  defined  in  Appendix  C: 
a.   General  formulas,  i  =  2,...,N: 

(A!)2 

0(1,1-1)   =   (u.  -u.  ,)•[ 1 7]  (1-1) 

1    1_1    48-(A.) 

l 

(A!)2      (A!) (A!  ,) 

Q(i,i)   =   (u.  -u.   )•[ 1 +  __i K^-] 

24-  (Ai)     16-  (Ai)Z 

(A!)2        3(A'   )(A!)     5(A'   )2 

+  (u.  .,  -u.)  •  [ ± j  +  i±i T-   +   ^-rl  (1-2) 

12.(Ai+1)     16.(A.+1)      48(A.+1)2 
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(A!)2        3(AJ. J  (A!) 
Q(i,i+1)   =   (u,^  -u,)-[ y  + 


1+1    X    12*^i+l)2     16*^i+l)2 

5(A!)2         _     _         (A •   )2 
-]  +  (u.^  -u.^)-[- 


2J     v  ui+?    i+l'   L  2 

48'Ui+l)  24*(Ai  +  l) 

(A!   )-(A!) 
+  — iii i_]  (1-3) 

(Ai+l)2 

0(1,1+2)   =   (u.,9-u.   )•[ — 5-]  (1-4) 

1+^    1+1    48-(Ai+1)^ 


b.   At  level  i  =  1 


5-  (u?  -  u, ) 
0(1,1)   =   |g (I"5) 


(A,)         _    _      (A!)2 

0(1,2)   =   (u9-u,)-[ 5-]  +  (u,  -u,)-[ j]  (1-6) 

1    12-  (A2)  24-  (A2) 

(AA)2 

0(1,3)   =   (u.  -u7)-[ 2]  (I_7) 

48- (A2) 


c.   At  level  i  =  N+l: 

(A'   )2 

Q(N+1,N)   =   (u    -u)-[ 5LJ: -.]  (1-8) 

48.(AN+1)2 
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—     —  2 

13-  (u.T,0  -uM,,  )  (A    ) 

Q(N+1,N+i>    =   —  48  +  (vrV-i — ^-^ 

48  ^+1    N    24-(AN+1)2 

+  — 5±i ^f  ]  (1-9) 

16-(AN+1) 


5.   Matrix  T: 

a.   General  formulas  for  i  =  3,...,N 

u.       -(A!)2  5-(A!)2  (A!)  (A!     ,) 

— ,  /   •        «       t   \  XX  X  r  X  X  XX 

T(l,l-1)        =       Tan-TK +    U-  '  tioo.  a + 


192-A.  i    l192-A.  24-A. 

l  ii 


(A!)2  A'            A!  u.       -(A! 

+  1 +       1+1 !_i   +  Jt±i ^_  (i-io) 

192-A.     .  192'A.     .       J  64-A.     .                            ^X    ±u; 

l+l  l+l  l+l 


u.     .  •  (A.)2  (A!)2  (A!)  (A!     .) 

m    ,     .  .     .  1-1  1  —  r  1  1  1-1 

T(l,l)  =  -T?To— A +      U/[C/!.A  + 


192-A.  i     L64-A.  48-A. 

l  ii 


h-(a:)2       7(a1+1)(a;)       17-(ai+1)2i 

192-A.+1  48-A.+1  192-A.+1 


17- (A!)2         11- (A'       )2         7- (A'       ) (A!) 

—  r  1  1+1  1+1 

+     ^  •      ,   '   [ttto-a +     inn.  a + 


i+1     L192-A.     .  192-A.     .  48-A.J_1 

l+l  l+l  l+l 

(A'        )2  (A'.     .)(A!)  u.     0-(A'       )2 

64-A.     .  48-A._Ll       J  192-A.     .  k 

l+l  l+l  l+l 
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T(i,i+1)   = 


ui-«Ai+l) 


64-A 


+  u.  ,,[ 


^i+l1 


i+1 


i+lL32-A 


(Aiti»(fli> 


i+1 


12-A 


i+1 


u    •  (A* 
i+2  l  i+l 


64-A 


i+1 


(1-12) 


b.   At  level  i  =  1: 


T(l,l)   = 


T(l,2)   = 


17-U.-A'    23-u2.A'    u3.A' 


192 


96 


192 
2 


U,-A'    13-u2*A^    U3(A3> 


64 


192 


64-A. 


(1-13) 


c.   At  level  i  =  2: 


T(2,'l)   = 


VA2    -   r13*A2    (A2)2     (A3)(A2)1 

+  u  •  f + +  : 

64      2  L  192     192-A,     24. A,   J 


u3-(A2) 


64-A. 


(1-14) 


T(2,2)   = 


ux-A'    _    7.A£    11- (A2)     7-(A3)(A2) 
+  u0-  [ + . —   + 


192 


2  L  192    192- A. 


48*  A. 


17-(A3)2     _    17-(A2)2    7,(A')2 

+  -z-z-x. : ]   +  U  •  [ : —   + 


192-A. 


3  L  192-A. 


96- A. 


(A'  )  (Al)     u,-(A!)2 
+  -^ — ]  + 


6.  A. 


192-A. 


(1-15) 


T(2,3)   = 


U2'(A3)2    -   r^3i2_    ^"ty.         U4'(A3)2 
64-A,    +  U3*L32-A,  +   12-A,   J  +   64-A, 


(1-16) 
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d.      At   level    i   =  N+l: 

um-(a'    )2  5*(\T.ui)2     (a;u  j(a;t) 

m/^.T     *t\                   N         N+l           ,     77             r             N+l  ,  N+l         N 

T(N+1,N)        =      -77-7 +    UATJ_n  •  [1Q-     , + 


64-AN+1  N+l     l192.AN+1  24.AN+1 

+    13'(AN+1}3    +    V2'^l 

192         J  64  v       x    ' 


u.-(A'     ,)2                         (A'    ,)2  (A'     ,)(A') 

m/%T,i    ».i\                   N         N+l           .     —             r       N+l  ,  N+l  N 

T  (N+l, N+l)      =      1qo.  A +   u..^,  •  [^^r^ + 


192-AN+1  N+l     L64.AN+1            48-A^ 

13  •  A  *  17 'U         • A ' 

,             /AN+1]  X/    UN+2    Vn        .                                            f,    18. 

+   —64 ]    +   192 U    18) 
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APPENDIX  J 
FORCED  VERTICAL  VELOCITY 

1.   The  contribution  to  the  surface  geopotential  from  the 
forced  vertical  velocity  is  $  , 


<j>M  sin  (^)sin(yx) 


t  <  T 


*s(x,t)   =   /  ,  (J-l) 

<|>M  sin(px)  t  >  T 


2   2 
where  <j>M  is  mountain  geopotential  (m  /s  ) ,  t  is  time,  and 

T  is  the  total  time  to  build  the  mountain.   d>._  is  a  constant, 

M 


*M   =   gHM  '  (J"2) 


where  g  is  gravity  and  HM  is  the  height  of  the  mountain. 
H„  is  a  parameter  specified  in  each  model.   H   =  750  meters 
in  the  thesis  experiments. 

2.   The  time  rate  of  change  of  $   is  separated  into  sine  and 
cosine  components  for  use  in  the  surface  geopotential  forecast 
equations, 


d$ 

3~   =   MTS1(t)  -cos(yx)  +  MTS2  (t)  -sin(yx)  ,  (J-3) 
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where  tt  =  tt  +  u  -  "5— .   Equation  (J-l)  is  substituted  into 
at    dt    sfc  dx 

equation  (J-3)  and  the  resultant  expression  is  separated 
into  sine  and  cosine  equations.   The  equations  to  calculate 
the  terms  MTS-l  and  MTS2  are 


and 


CS1(t)   =  ,  (J-4) 

sfc  YM 


— - — sm(^)-cos(^)    t<T 
CS2(t)   =  <  .  (J-5) 

0  t  >  T 


These  terms  are  calculated  for  each  time  step  in  the  model's 
forecast  subroutine. 
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APPENDIX  K 
DIABATIC  HEATING  TERMS 

For  the  diabatic  heating  function  defined  in  equation 
(3-1),  the  maximum  heating  occurs  at  Z  =  0.50,  the  midpoint 
of  the  heated  layer.   Temperature  and  diabatic  heating  are 
defined  for  grid  A  at  the  staggered  levels,  while  they  are 
defined  at  the  unstaggered  levels  for  grids  B  and  C.   Conse- 
quently, the  rate  of  heating  differs  between  the  staggered 
and  unstaggered  models.   In  these  experiments,  the  heated 
layer  is  between  Z  =  0.40  and  Z  =  0.60,  the  heating  rate  is 
5.0  K/day,  and  only  the  cosine  term,  Ql,  is  nonzero  in  the 
heated  layer.   The  value  of  the  heating  term  is  listed  below 
for  staggered  and  unstaggered  levels  for  six,  twelve  and 
sixty-layer  models. 

Grid  A  60-Layer  Models      Grid  B  and  C  60-Layer  Models 

Z_  Ql 

0.400  0.0 

0.417  0.388E-05 

0.433  0.145E-04 

0.450  0.289E-04 

0.467  0.434E-04 

0.483  0.540E-04 

0.500  0.579E-04 

0.517  0.540E-04 

0.533  0.434E-04 

0.550  0.289E-04 
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z 

Ql 

0.392 

0.0 

0.408 

0.986E-06 

0.425 

0.847E-05 

0.442 

0.214E-04 

0.458 

0.364E-04 

0.475 

0.493E-04 

0.492 

0.569E-04 

0.508 

0.569E-04 

0.525 

0.493E-04 

0.542 

0.364E-04 

0.558  0.214E-04 

0.575  0.847E-05 

0.592  0.986E-06 
0.608         0.0 


0.567 

0.145E-04 

0.583 

0.388E-05 

0.600 

0.0 

Grid  A  12-Layer  Models 


Grid  B  and  C  12-Layer  Models 


z 

91 

Z  . 

Si 

0.375 

0.0 

0.333 

0.0 

0.458 

0.364E-04 

0.417 

0.388E-05 

0.542 

0.364E-04 

0.500 

0.579E-04 

0.625 

0.0 

0.583 

0.388E-05 

0.667 

0.0 

Grid  A 

6-Layer  Models 

Grid  B 

and  C  6-Layer  Models 

Z 

Qi 

Z 

Qi 

0.250 

0.0 

0.333 

0.0 

0.417 

0.388E-05 

0.500 

0.579E-04 

0.583 

0.388E-05 

0.667 

0.0 

0.750 

0.0 
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